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Chapter 1 


Introduction 


1.1 Outline 

Control of air contaminants is a crucial factor in the safety considerations 
of crewed space flight. Indoor air quality needs to be closely monitored during long 
range missions such as a Mars mission, and also on large complex space structures 
such as the International Space Station. This work mainly pertains to the detection 
and simulation of air contaminants in the space station, though much of the work is 
easily extended to buildings, and issues of ventilation systems. 

Here we propose a method with which to track the presence of contaminants 
using an accurate physical model, and also develop a robust procedure that would 
raise alarms when certain tolerance levels are exceeded. A part of this research 
concerns the modeling of air flow inside a spacecraft, and the consequent dispersal 
pattern of contaminants. Our objective is to also monitor the contaminants on-line, 
so we develop a state estimation procedure that makes use of the measurements 
from a sensor system and determines an optimal estimate of the contamination in 
the system as a function of time and space. The real-time optimal estimates in turn 
are used to detect faults in the system and also offer diagnoses as to their sources. 

This work is concerned with the monitoring of air contaminants aboard 
future generation spacecraft and seeks to satisfy NASA’s requirements as outlined 
in their Strategic Plan document (Technology Development Requirements, 1996). 
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Most of this work was undertaken to satisfy the requirements of NASA’s Advanced 
Environmental Monitoring and Control Program, with a view to developing an in- 
telligent monitoring system for Space Station missions. The current NASA Strategic 
Plan has as one of its stated goals “to explore, use, and enable the development of 
space for human enterprise”. The goal is to be accomplished in three time periods 

• 1996-2002: Establish permanent presence in low Earth orbit by constructing 
and using the ISS, 

• 2003-2009: Operate the ISS cost effectively, with a subgoal to “achieve ad- 
vanced life support systems to close spacecraft air/water loops,” and, 

• 2010-2020 and beyond: Conduct international human missions to planetary 
bodies in our solar system. 

Even though this work is targeted at future generation spacecraft and space sta- 
tions, many of the specifications used in this work pertain to the International Space 
Station, scheduled for launch in late 1998. The International Space Station is a col- 
laborative effort with participation by the governments of the United States, Canada, 
Europe, Japan, and Russia. The configuration will include a Hab and a Lab Element, 
two Nodes, and two International modules (the European Space Agency Attached 
Pressurized Module and the Japanese Experimental Module). Other relevant mis- 
sions where this work applies include the manned missions to Mars, the Mars Short 
Visit, the Mars Human -Tended Outpost, and ;he Mars Permanent Outpost, where 
astronauts are expected to spend up to 600 dajs, and where the luxury of returning 
to earth for a cleanup in the case of a contaminant leak does not exist. Contaminants 
that are to be monitored include carbon dioxide, carbon monoxide, and volatile or- 
ganics. According to NASA, primary chemicals of concern are nitrogen tetroxide, 
monomethyl hydrazine, ammonia, and Halon 1301. Studies aboard the Mir Sta- 
tion(Cole et al., 1996) have shown that about 45 compounds (32 target compounds 
and 13 non-target compounds) were consistently detected in air samples that were 
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collected during the missions, though none of the compounds were present at toxic 
levels. Still, contaminant monitoring remains an important focus area in ensuring 
the safety of humans in space. 

Major sources of contaminants in the space station include off-gassing of 
cabin materials and hardware, use of utility chemicals, and metabolic waste products 
of crew members. Minor sources of contaminants include electrical equipment, mi- 
crobial metabolism, leakage during experiments using chemicals, leakage from envi- 
ronmental or flight control systems, volatile food components, and reaction products 
from the environmental control and life support systems. Table 1.1 lists some of the 
substances being monitored and their sources aboard the spacecraft. 

The National Research Council (National Research Council, 1984), in its 
various studies, has prescribed (National Academy of Sciences, 1981) spacecraft max- 
imum allowable concentrations (SMACs), which are not to be exceeded under any 
circumstances (National Research Council, 1996). These concentrations are based 
on studies that link contaminant concentrations to the impairment of normal human 
activities. 

The fault detection and diagnosis system is a synthesis of different math- 
ematical procedures, which are functionally independent, but which come together 
to provide an effective structure that serves the purpose of monitoring the presence 
of airborne contaminants. 

1.2 Space Station Environmental Control and Life Support System 

In this section, the general layout of the International Space Station, and 
some of the components of the Air revitalization system are described. 

Future missions, especially the long range missions, will increasingly have 
to be materially closed systems, since the cost of carrying spare oxygen would be 
prohibitive. Future missions will also involve growing food, and the complexity of the 
revitalization system will increase to account for many more possible contaminants. 
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Table 1.1: Some commonly monitored substances and their sources aboard spacecraft 
(National Research Council, 1992) 


Monitored substance 

Spacecraft source 

Oxygen 

Required component in cabin air 

Carbon dioxide 

Required component in cabin air 

Carbon monoxide 

Product of incomplete combustion 

Nitrogen 

Inert component in air 

Halon 

Diffusion from the Shuttle to the Station 

1-Butanol 

Off-gassing from flight hardware 
and from human metabolism 

Diacetone alcohol 

Off-gassing from paint that is not totally cured 
and from hardware off-gassing 

Ethanol 

Cleaner /disinfectant use 

Ethyl benzene 

Off-gassing from nonmetallic materials 

Ethylene Glycol 

antifreeze and off-gassing 

Glutar aldehyde 

Cell-Tissue fixatives 

Trichloroethylene 

Off-gassing 

Xylene 

Off-gassing 
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In addition, NASA’s Advanced technology requirements outlines certain features 
that must be included in future life support systems. Some of those features that 
this work attempts to demonstrate include: 

• Ambient air in the cabin must be monitored at selected locations, every 15 
seconds, for the species O 2 , CO 2 , and CO. 

• Toxicity of air must be reported in terms of specific major and trace species 
concentrations and their rates of rise. 

• A computer model shall be available to predict the behavior and the contami- 
nation removal capabilities for contaminants that could suddenly be released 
into the atmosphere. The model must be experimentally verified and be ca- 
pable of spatial resolution to the module level and temporal resolution to 
0.5 hour. 

• major air components shall be monitored on a near-continuous basis in the 
habitat atmosphere 

The Environmental Control and Life Support system is divided into seven 
major subsystems, the temperature and humidity control (THC), atmosphere control 
and supply (ACS), atmosphere revitalization (AR), fire detection and suppression 
(FRS), water recovery and management (WRM), waste management (WM), and the 
Vacuum System (VS). 

Their functions include (Reuter, 1998) : 

• Atmosphere Revitalization: 

* Control and disposal of carbon dioxide 

* Control of airborne trace contaminants 

* Oxygen (O 2 ) supply via generation 
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* Atmospheric monitoring of primary constituents, including 0 2 , C0 2 , 
nitrogen (N 2 ), hydrogen (H 2 ), methane (CH 4 ), and water vapor 

* Airborne particulate and microbial control 

• Temperature and Humidity Control 

* Cabin air temperature and humidity control 

* Equipment air cooling 

* Inter- and intra-module ventilation for crew comfort and station level 
control of 0 2 , C0 2 , and trace contaminants 

• Atmosphere Control and Supply 

* Total pressure and 0 2 partial pressure control 

* Total pressure monitoring and monitoring of loss of pressure 

* Stored gaseous N 2 and 0 2 supply and replenishment 

• Fire Detection and Suppression 

* Smoke detection 

* Fire extinguishment 

• Waste Recovery and Management 

* Potable and hygienic water supply 

* Waste water and urine water collection, recovery, and disposal 

• Waste Management: 

* Urine/fecal collection and recovery 

• Vacuum System 


* Vacuum venting and maintenance for payload support 
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Some of the air and water quality requirements which are to be maintained 
in the ISS are shown in Table 1.2. 

A test of the Life Support System was carried out in the Johnson Space 
Center’s Integrated Life Support System Test facility, in which a four-member crew 
spent thirty continuous days in a closed chamber. All test objectives were accom- 
plished (Lunar-Mars, 1997), and no SMACs were violated. 

1.3 Previous work 

Space environment monitoring has been in place for as long as there have 
been crewed missions, with the levels of sophistication changing with the complexity 
and duration of the missions, and along with the developments in computational 
and sensor technologies. Two computer models that represent the present genera- 
tion of Space Environment monitoring models are the Trace Contaminant Control 
Simulation (TCCS) (Perry, 1993) and the Computer Aided System Engineering and 
Analysis (CASE/A) (CASE/ A, 1990) modeling package. Both the packages model 
the space station modules as well-stirred tanks, where each module is represented by 
its average concentration. The CASE/ A package is more flexible and user-friendly, 
and provides a means for simulation of a number of interconnected well-stirred tanks. 
Both the packages suffer from the limitation that stagnation points within the cabin, 
and non-uniform forced convection patterns cannot be represented within the model. 
A study of the inhalation risks aboard spacecraft (Todd et al., 1994), where the 
Space Station was modeled as a series of well stirred chemical reactors, improved 
on the lumped analysis by providing more information about spatial variations of 
contaminants aboard spacecraft. A comprehensive study of the sensor system was 
undertaken (Smith, 1996), which used a lumped system of modeling the transport 
in order to optimize the sensor selection process. 

The first significant shift away from lumped analysis in contaminant dis- 
persal modeling was in the work of Skliar and Ramirez (Skliar and Ramirez, 1997a), 
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where each cabin module was represented using a two dimensional mesh and a solu- 
tion to the reaction-diffusion equation was obtained using a finite difference scheme. 
This approach has the advantage of providing information in two spatial directions. 

Blackwell (Blackwell, 1998) has proposed a a fault detection and location 
procedure based on physical laws and modeling, analysis, and simulation. His work 
is based on optimal control theory, and uses analytical redundancy for the detection 
part, but the work does not specifically address the physics of the problem, and 
merely lays out the structure for the procedure. 

1.4 Salient characteristics of this research 

This research seeks to build on the two-dimensional model developed by 
Skliar (1996), and extend that work to three dimensions. It also seeks to use a more 
rigorous computational fluid dynamics solution to the flow equations that what has 
been used hitherto for this purpose, one that vastly reduces the modeling error. A 
fault detection algorithm is implemented with the ability to distinguish between pro- 
cess and sensor faults. The fault detection algorithm uses the concept of analytical 
redundancy to detect the faults. 

The final part of the research focuses on diagnosing the fault in the system, 
primarily that of an unknown contaminant source in the cabin. This is an inverse 
problem that is ill-posed and has no unique solution, so the attempt here will be 
to obtain a satisfactory solution to a required degree of accuracy. To this end, an 
Extended Implicit Kalman Filter is developed, which is an extension of the Implicit 
Kalman Filter suitable for joint state and parameter estimation. The filter essentially 
maintains the same structure as the original filter, and many of the algorithms 
remain the same as before. The filter requires an initial guess for the location and 
capacity of the unknown source, for which purpose we use pre-calculated sensitivity 
matrices that contain information about the local variation in concentration for 
perturbations throughout the cabin. The combination of these two techniques makes 
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for an accurate and efficient algorithm for obtaining a solution to the inverse problem. 

The proposed procedure for air contaminant monitoring is shown in Fig. 
1.1. The same procedure can be applied for a variety of substances,, and some 
algorithms can even be used to monitor air pressure and some other parameters. 

1.5 Outline of Report 

This report is organized as follows. Chapter 2 discusses the flow modeling 
work that we performed in order to simulate the air flow aboard the International 
Space Station, and Chapter 3 describes the mathematical modeling of the contam- 
inant dispersal process along with its numerical solution, and sample profiles that 
were obtained. The State estimation procedure is discussed in Chapter 4, and the 
use of the State Estimation Procedure for fault detection is described in Chapter 5. 
Chapter 6 contains the discussion of the fault diagnosis algorithms that are proposed, 
and finally Chapter 7 presents the conclusions and the significant contributions of 
this research. 
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Figure 1.1: Air contaminant monitoring system 







Chapter 2 


Indoor air flow 


2.1 Literature review 

Air flows inside enclosures have been a subject of active interest within the 
building systems community. Computational fluid dynamics (CFD) has been used 
for predicting room air movement since the seventies. There even exists a public 
domain software (Kurabuchi et al., 1990; Said et al., 1995) called EXACT3, which 
is a three dimensional finite difference computer program for simulating buoyant 
turbulent airflow within buildings. In recent years, much effort has been made to 
enhance CFD as a reliable tool for the evaluation of air flows. CFD has been used 
in studying clean room air-flows (Kuehn, 19£8; Yamamoto et al., 1988) because 
of the need to keep the clean room free of particles and particulates, and air flow 
becomes a critical parameter in such cases. Space application in the context of 
contaminant dispersion is very similar in scope to clean room applications, though 
there is little work reported in the literature that pertains to space applications. 
Recently (Tam, 1998), an interesting study eva uated the application of CFD in the 
software design of environmental control and life support systems, and found that 
atmospheric conditions within the Space Statio l could be adequately modeled using 
the Resource Utilization Planning and System Modeling (RUPSM) scheme. 

Most of the work in the area has remained computational, though a few 
validating experimental results also exist in the literature. The International Energy 
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Agency, though a subsidiary research group, Annex 20, measured velocities, temper- 
atures, and turbulence velocity scales in full scale rooms (Whittle and Clancy, 1991). 
Canadian researchers (Barber et al., 1982) did studies on the correlation between the 
velocity of the inlet jet and the floor velocity, and proposed a jet momentum number 
that would measure the energy contained in the diffuser jet relative to the room air 
volume. In this work, we use CFD to provide us with the information regarding the 
flow which is then used as an input to the mathematical model for the diffusion and 
for the procedure that estimates the current concentration of contaminants in the 
cabin. This development marks an important step in our estimation procedure since 
the accuracy of the procedure is largely dependent on the accuracy of the flow model 
since most of the transport is occurring through convective diffusion. 

A detailed knowledge of the flow field is required in order to ensure that 
the ventilation system is performing adequately, and to provide information about 
local variations in the concentration profile of the contaminants. Another advantage 
in using CFD is that it enables the calculation of quantities like turbulence intensity 
which have direct effects on the comfort level of people inside the cabin. CFD is an 
inexpensive tool for such studies and has been used to study the effect of different 
ventilation techniques (Gan, 1995) on thermal comfort in buildings. Research has 
shown that lower turbulence intensities contribute to higher comfort levels (Zhang 
et al., 1992). 

2.2 Air Flow Modeling 

We assume that the air flow is steady and incompressible, and solve the 
three-dimensional Navier-Stokes equations for the appropriate boundary and initial 
conditions. Air under atmospheric pressure and for the low velocity flows that are 
characteristic of room air flows is expected to be incompressible, and we invoke the 
steady state assumption because solving for transient cabin flows is computationally 
too intensive to be used in a real-time operation. We tried simulating the air flow 
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for both 2-D and 3-D cases. Earlier work used a two-dimensional flow field (Skliar 
and Ramirez, 1997a) as an input to their state estimation procedure. In a significant 
development in this work, we consider modeling the three-dimensional geometry of 
the space station module. The 2-D case has the advantage of requiring far less 
computational time, whereas it suffers from a lack of information about the third 
dimension. 

Air flows inside enclosures are usually turbulent, random, and highly recir- 
culating (Zhang et al., 1992). In this work, we solve the equations for both laminar 
and turbulent cases. The geometry chosen here follows experimental test used in 
previous studies (Son and Barker, 1997) and accurately represents the US Space 
Station Lab module. We used this geometry for our simulations in order to have an 
experimental set of results to validate our simulations. The simulated cabin (Fig. 
2.1) is 6 m long, 2 m wide, and 2 m high (approximately. 20’ X T X 7’). There 
are two inlets and two outlets for the air. The Temperature and Humidity Control 
(THC) is the primary air supply which supplies regulated air into the cabin and 
is one of the primary subsystems in the Environmental Control and Life Support 
Systems for the Space Station (see Son and Barker, 1993). The Intermodule ven- 
tilation (IMV) airflow assemblies are used to interchange airflow between modules. 
One would expect that the THC air is relatively contaminant free since it is filtered, 
whereas the IMV could have trace contaminants generated in other modules, both 
routine contaminants and those released due to accidents. 

The air-flow model is based on the ccntinuity equation, the Navier-Stokes 
equation, thermal energy equation, and the concentration equation together with 
the k — e turbulence model equations, for the case of turbulent flow. The k — e model 
(Anderson et al., 1984; Whitaker, 1968) uses the kinetic energy of turbulence /c, and 
its dissipation rate e as the two scales. This introduces two additional transport 
equations. In the k — e model, the turbulent viscosity fi t also known as the eddy 
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Figure 2.1: Sketch of the model cabin. 




16 


viscosity is defined by the relation: 


Pt = c^pk 2 ie (2.1) 

where is an empirical constant, and p is the density of the fluid. 

For a rectangular cabin geometry, the equations of continuity, momentum, 
energy and mass for an incompressible flow are as follows. The symbols in the 
equations are defined in the nomenclature section at the end of this thesis. 
Continuity: 


V • u = 0 


( 2 . 2 ) 


Momentum: 


The momentum equation for the three co-ordinate directions is given below. 


x-direction 


d(pu) d(pu 2 ) d(puv) d(puw) 

dt dx dy + dz 


(2.3) 


dp d _ du. d . ,dv du., d . ,du dw,, 

■fe + te (w ' ° + ^ + ay* W 8* + + fo 1 + a* )] + pU 


y-direction 


d{pv) d(puv) d(pv 2 ) d(pvw) 

dt dx dy dz 


(2.4) 


dp d . .dv 

dy dx ^ dz 


. ^ u \i , & /\x-f r, & v \ d . dw dv.. , 

+ a ? 11 + aj? ( 2 "a» ) + af + +p/ » 


z-direction 


9(/nu) c?(/cmiu) d(pvw) 

r-v ^ 




dx 


+ 


d(pw 2 ) 


’ <9z 

dp d x ( du dw d r ( dw dv . d 3u\ 

“S + & w ai + fe ’ 1 + + fc 11 + fe (AV ■ u + 2p fe) +p/ 


(2.5) 
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Energy: 

The energy equation is given below for the sake of completeness. The 
Space Station environment is expected to be isothermal for the most part, and 
minor temperatures would not affect the flow field much since buoyancy is not an 
important factor under micro-gravity. 


D. u 2 . . d n dT. d n dT. d dT 

p Wt { ‘ + Y } = pq + di (k te> + &<%> + &<*& > 

.du dv dw , w <9u dv dw.n 


(2.6) 


# <£> 2 + + ^) 2 + (£ + + & + £) 2 + (= + S 2 ! 


<9y' 




k (9x ' 


, dv 
' dz 


du k 
dy 


2.3 Solving the Navier-Stokes equations 


The equations were solved in a coupled manner using the Fluid Dynam- 
ics Analysis Package (FIDAP Version 7.62)(FIDAP, 1993). A finite element mesh 
grid was developed for the two-dimensional and the three-dimensional problems with 
specified nodal boundary conditions. An eight-node brick was used as the basic finite 
element for the purpose of discretization. The velocity components were approxi- 
mated using trilinear interpolation functions within the elements. The pressure was 
discretized in a piecewise continuous manner, with the pressure degree of freedom 
associated with the element centroid. A segregated solver was used to solve the re- 
sulting non-linear equations. The segregated solver decouples the equations for the 
purpose of solution, and sequentially solves them, using the results of one equation 
in the next, and so on. This increases the CPU time needed, but conserves memory, 
and has been found to be very useful for large mesh sizes. A variety of boundary con- 
ditions was tried, though for the sake of conciseness, only two cases will be discussed 
here. The Convergence criterion required the residuals of the equations (velocity 
components and the pressure) to be below 0.0001. Most of the computations were 
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performed on the SGI Power Challenge Array, in a parallel mode with either two 
or three processors, available to us through the National Center for Supercomput- 
ing Applications, Urbana-Champaign. The runs required varying times to converge, 
approximately in the range of 60 CPU hours for each. 

2.4 Flow field results 

Case 1: Laminar Flow- high flow rate 

The geometry has already been described. For this case, we used an IMV 
with a duct velocity of 4.8 m/s and a THC with a velocity of 6.0 m/s. The cabin is 
assumed to be isothermal. The ducts leading to the outlet are modeled using free 
boundary conditions, i.e. the values of the velocities are allowed to float to satisfy 
the Navier-Stokes equations. The no slip boundary condition was invoked at all the 
walls. We assume that the flow is laminar, and that it is a steady flow field. 

Figures 2.2 - 2.4 illustrate the flow profiles that we obtained. The contour 
of the magnitude of the velocities and the velocity vectors themselves are depicted. 
We show three horizontal slices of the box, one near the top and bottom, and one 
halfway up the cabin to illustrate the variations. The slice near the bottom is mainly 
dominated by the exit of the THC duct. Note that the flow leaves at an angle to 
the duct due to the blast of air that blows in the x-direction. The slice near the 
center clearly shows the profile near the inlet and outlet for the IMV ducts. The 
flow spreads out throughout the room. The laminar case shows no recirculation cells. 
The slice from the top of the cabin shows the low entering the cabin, and the cells 
formed as the jet curves downward. 

Case 2: Turbulent Flow- high flow ral e 

Turbulence, in a sense, is still an unsolved problem. The presence of a 
length and a time scale much smaller than the physical problem presents a scenario 
where the exact solution to equations cannot be obtained. In addition, for room and 
cabin flows, it is difficult to predict the onset of turbulence. A statistical approach 
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Figure 2.3: Speed Contour and Velocity Vector for the middle slice under laminar 
flow conditions for a THC flow of 6.0 ms -1 and an IMV flow of 4.8ms -1 . The speed 
and velocity are expressed in the units of ms -1 . 
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Figure 2.4: Speed Contour and Velocity Vector for the top slice under laminar flow 
conditions for a THC flow of 6.0 ms -1 and an IMV flow of 4.8ms” 1 . The speed and 
velocity are expressed in the units of ms -1 . 
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is usually used, and the equations are averaged over a time scale that is long com- 
pared to that of turbulent motion. The resulting averaged equations then describe 
the distribution of the mean velocity, pressure, temperature and the other quanti- 
ties of interest. Detailed derivations of the equations can be found in any advanced 
book on fluid motion (Anderson et ah, 1984). We use the two-equation model, 
briefly touched upon earlier. For isothermal flow with no mass transfer, the recom- 
mended set of model parameters were used for these empirical constants defined in 
the nomenclature section. 

= 0.09, Ok = 1.00, <x c = 1.30, ci = 1.44, c 2 = 1.92 (2.7) 

Figures (2. 5-2. 7) are graphical representations of the turbulent flow simulations. 

Case 3: Laminar Flow -Low flow rate 

For the same geometry as before, the flow rates were decreased, and the 
steady state flow profiles were recalculated for an IMV flow rate of 0.15 m/s and a 
THC flow rate of 0.3m/s. All other conditions were maintained at previous levels. 
The profiles are shown in Fig. 2.8-2.10. 

Case 4: Turbulent Flow -Low flow rate 

The low velocity flow field calculation was repeated for turbulent flow. Fig- 
ures 2.11-2.13 show the contours for the profiles that were obtained. 

No major differences were noticed in the flow profiles obtained for high and 
low velocity duct flows. The patterns of flow essentially remained the same. A more 
thorough study of inlet velocities and their effect on cabin air flows is in order but 
beyond the scope of this work. The turbulen flow profiles closely resemble those 
obtained previously in experiments (Son and Barker, 1997) in the Space Station 
simulator. The differences can be attributed to the minor differences in the geometry 
in the region of the hatches connecting the modules, and the roundedness of their 
hatches. The existence of recirculation cells is the significant difference between the 
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Figure 2.5: Speed Contour and Velocity Vector for the bottom slice under turbulent 
flow conditions for a THC flow of 6.0 ms" 1 and an IMV flow of 4.8ms -1 . The speed 
and velocity are expressed in the units of ms -1 . 
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Figure 2.7: Speed Contour and Velocity Vector for the top slice under turbulent flow 
conditions for a THC flow of 6.0 ms -1 and an IMV flow of 4.8ms -1 . The speed and 
velocity are expressed in the units of ms -1 . 
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Figure 2.11: Speed Contours for the top slice - (1.95m from the floor) under turbulent flow conditions for a THC flow of 0.3 
ms -1 and an IMV flow of 0.15ms" 1 . The speed and velocity are expressed in the units of ms” 1 . 








Figure 2.12: Speed Contours for the middle size (0.9m from the floor) turbulent flow conditions for a THC flow of 0.3 
an IMV flow of 0.15ms *. The speed and velocity are expressed in the units of ms -1 . 
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Table 2.1: Normal mesh-errors for the velocity components and the pressure equa- 
tions. 


Degree of Freedom 

RESIDUE 

u 

0.72617E-03 

V 

0.10399E-02 

w 

0.21082E-03 

P 

0.39848E-03 


laminar and the turbulent profiles. We believe that the actual flow in the Space 
Station is mostly turbulent, and the the profiles we obtain are characteristic of low- 
velocity turbulence flows. 

2.5 Mesh refinement studies 

Numerical simulations are of course, subject to errors, and are closely re- 
lated to the coarseness of the mesh used in the simulations. One common way of 
validating CFD simulation results is to refine the mesh being used and noting that 
there was no major change in the solution obtained for the same geometry, initial 
and boundary conditions. Tables 2.1 and 2.2 show the error residuals for each of the 
four equations, representing the velocities in the three co-ordinate directions and the 
pressure. Results are shown for both a normal mesh and for a refined mesh (dou- 
ble the number of mesh points). The results indicate that further refinement will 
nor substantially change the overall flow profile. The flow profiles obtained in the 
previous sections suggest that a refinement of the mesh near the flow inlets would 
improve the accuracy of the solution near the inlets. 

2.6 Complex flows 

The flow in a cabin under operation is likely to be different from what we 
obtain by simulating an empty cabin, as we have done until now. Humans in the 
cabin and equipment are likely to cause more turbulence and recirculation cells. A 
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Table 2.2: Refined mesh-errors for the velocity components and the pressure equa- 
tions. 


Degree of Freedom 

RESIDUE 

u 

0.12671E-02 

v 

0.14342E-02 

w 

0.10849E-02 

p 

0.71476E-03 


detailed analysis of complex flows is out of scope of this work, but we attempt to 
obtain some preliminary results on what flow fields look like with objects in the 
cabin. 

A set of simulations were carried out for the geometry shown in Fig. 2.14. 
Speed contours for three horizontal slices are shown in Figs. 2.15-2.17. The turbu- 
lence in the cabin was found to increase with the presence of an object, with more 
dead zones. The number of recirculation cells has increased. The velocities did not 
vary too much, due to the low velocity of the air flow. 

2.7 Summary and Conclusions 

Our objective in studying the air flow inside the cabin was to arrive at a 
basic understanding of cabin flows aboard the Space Station, and to obtain a few 
sample air flows that could be used as an input to our transport model, developed 
in the next Chapter. The flow profiles are very important to the transport model 
since most of the mass transfer in the cabin occurs as convective transport, and 
the accuracy of the flow field will therefore largely control the accuracy of the final 
transport model. Space Station flows are very poorly understood at present, and 
considerable further work is needed, both simulation of flows using CFD techniques, 
and experimental work which can validate those simulations. 
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Figure 2.14: The cabin with an object inside, showing the THC inlet flow of 0.3 ms and the IMV outlet flow. The velocity is 
expressed in the units of ms -1 . 




Figure 2.15: Speed Contours for the top slice - (1.95m from the floor) with an object (not shown) in the cabin under turbulent 
flow conditions for a THC flow of 0.3 ms' 1 and an IMV flow of 0.5ms -1 . The speed and velocity are expressed in the units of 
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Chapter 3 


Air contaminant dispersion modeling 


3.1 Mathematical model 

One of the goals of this work was the development of a complete three- 
dimensional transport model for the air-borne contaminants. The essential structure 
of the model is shown in Figure 3.1. Previous work used either the control- volume 
approach, in which the cabin was modeled as a collection of well-mixed reactors 
(National Academy of Sciences, 1981), or a space-averaged two-dimensional model 
(Skliax and Ramirez, 1997a). 

Nazaroff reports a study of the effects of indoor air pollutants in which the 
indoor air was tagged with a non-reactive tracer and the decay of its concentration 
was monitored. They also extended this work to the study of flows between rooms 
using a combination of remote sensing and computed tomography techniques, which 
yields accurate results for the dispersion of air pollutants. 

In this section, we present our model o'the transport process, the discretiza- 
tion scheme to convert the partial differential equations to discrete representations, 
and discuss our solution scheme. The three-dimensional transport model developed 
here essentially extends the two dimensional model of Skliar (1996). 

Assuming a constant density of spaceci aft atmosphere and a constant molec- 
ular diffusivity, Dm, the differential mass balance of the air-borne contaminant 
(Bird et al., 1960) with concentration q results in the following fundamental three- 
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Figure 3.1: Inputs and output from the air contaminant transport model 
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dimensional convection-diffusion transport model 

do 9 

t^ + u- Vq = D M V 2 q + f (3.1) 

where u is the bulk air velocity, and / is the function that describes the capacity 
and location of contaminant sources and sinks. In this work, we are assuming that 
the contaminant is passive, i.e. it is transported with air at the same velocity in the 
field. In addition, we are assuming that the contaminant undergoes no chemical or 
physical transformations during its transport. 

Equation 3.1 is applicable to both laminar and turbulent flow. However, 
in the case of turbulent flow, the velocity vector is extremely random, and so we 
resort to using the time averaged equations instead. The idea is to average the 
Fickian model over a time interval long enough for the integral of the instantaneous 
fluctuations to become zero. For the case of turbulent flow, therefore, we treat both 
the flow velocity and the concentration q as stochastic quantities. The transport 
equation, for the case of turbulent flow is written as 

+ U ■ Vg = S7 2 D M q + / (3.2) 

where the overbar indicates that they are time -averaged quantities. 

The eddy diffusivity, Dm, which is tae diffusivity under turbulent condi- 
tions, is a function of the flow field and is therefore not uniform throughout the 
geometry. 

3.2 Computer Implementation of the three Dimensional Model 

We solve the model Eqs. 3.1 and 3.2 above using a simple finite differenc- 
ing scheme. In this section, we discuss how we discretize the equations, and then 
outline the solution technique used to solve the same equations. Since the flow field 
is expected to be turbulent in most cases, this derivation will proceed under that 
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assumption, and using the eddy diffusivity as the diffusivity coefficient. If laminar 
conditions are found to prevail, the molecular diffusivity should be used, which is 
independent of the flow field, and therefore is uniform throughout the cabin. 

The diffusive terms are discretized using a second center difference scheme 
and appear as 


d TT dq - A 
Yx M dx ~ dn+hp ’ r 


<?n+l,p,r Qn y p y r A Qn y p y r Qn~l,p,r 
' ” a n,p y r - rt 


Ax 2 


Ax 2 


( 3 . 3 ) 


d dq 
dy dy 


. gn,p+l,r Qn,p,r _ , 

a n,P+ l.r ^y2 


x n,p,r 


Qn,p } r Qn,p-l,r 

Ay 2 


dz dz 


, Qn,p,r+\ Qn,p,r , 9n,p,r 9n,p,r- 1 

“n,p,r+l Az 2 a n,p,r 


( 3 . 4 ) 

( 3 . 5 ) 


where Ax, Ay ,A z are discretization steps along coordinates x,y, and z, d is the 
discrete analog of the diffusivity and the subscript is used to specify a point on the 
spatial mesh {(n,p,r) | n = 1, N,p = 1 ,P,r = 1 ,#}. N,P, and R determine the 
mesh size used in discretizing the cabin geometry. 

The convective terms are discretized using the upwind differencing scheme 
in order to eliminate any possible oscillatory effects in the solution. The convec- 
tive terms for the East-West, South-North, and Up-Down directions are as follows. 
East- West 


where 


udq 




dx 

^n,p,r 

Ax 

•> 

, e _ J 

\ Qn,p,r 

if ^n,p,r 

V 

p 

n,P,r j 

l Qn-l y p y 

r if ^n,p, 

r < 0, 

W _ J 

[ 9n-l,p, 

r if u nj p 5 

r > 0, 

n ,p,r j 
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if w n y p,r 

< 0. 


( 3 . 6 ) 


South-North 


o S N 

VOq 9n,p,r Qn y p y r 

dy ~ Vn ’ p ' r Ay ’ 


( 3 . 7 ) 



42 


where 


Up- Down 


where 


c f 9n,p-l,r if u 7i,p,r > t 
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l 9 n,p,r if ^a,p,r < 0 , 
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?n >P> r ] . r 
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= I 


(3.8) 


Qn,p,r— 1 if w n,p,r > 0, 

Qn,p,r if w i,p,r < 0. 

The application of the center difference approximation of the time derivative yields 
the following discrete analog of the three dimensional transport model. 


0 m+l _ a m 

^n,p,r **n,p,r 

~ A t 


u 


7i,p,r 


V 


Ax Un ' p ’ r 


Ax 


, r [4 r r +1 -[C] mfl , Kr-fcr] 

+ n ' p ' r A y +Vn *’ r A y 

[g^P.r] ~ [gn, p,r] [^p,r] ~ [tfn,; 

^ ^n^p.r — 


PP 


Az 


Az 


r7 m + 1 — i 7 m +i „rn+l _ n m + 1 

j <*n+l,p,r Vn,p,r j *in,p,r Vn-1 ,p,r 

dn+1 ’ p ' r A^ ^ A^ 


rt 771 , __ n m n m _ „m 

, ^ Vn+l,p,r Vn,p,r j ^/n,p,r 9n-l,p,r 

“T a n+l,p,r A _2 a np,r 


Ax 2 "" p ’ r 

a m + 1 _ „m+l m+1 _ _m+l 

j “n,p+l,r Vn.p.r j v n ,p,r yn,p-l,r 

n ’ p+1 ’ r Ay 2 n>p>r Ay 2 


< 7 m , — n m n m — n m 

■ ^ "n f p+l,r Vn,p,r j *n,p,r Vn,p-l,r 

+ ‘t.J+l.r ^5 <*»«• ^2 


m+1 __ m+1 „m+ 1 _ „m + l 

W ^n^p,r+l Vn,p,r ^ *n,p,r Vn,p,r-1 

a n,p,r+l 2 «n,p,r ^2 



43 


+ 


ri^ n^' 

j H n,p,r + 1 'm.p.r j ^n.p.r Vn.p, r-1 , xm +l 

«n,p,r+l ^2 “”>P. r ^ z 2 ^ Jn,p,r 


(3.9) 


where represents the value of the time dependent source function / evaluated 
at the current time step, At is the time discretization step, and the superscript 
m = 0, 1, 2, . . . is used to identify an instance t = (m + 1) At for which the solution 
of the equation is sought. 


3.2.1 Numerical solution of the transport model equations 

Because of their poor stability properties, explicit difference methods are 
rarely used to solve initial and boundary- value problems in two or more space dimen- 
sions. The solution scheme used here is the classic Alternating Direction Implicit 
scheme (ADI) (Douglas and Rachford, 1956; Douglas, 1962), which invokes the prop- 
erty of operator splitting and converts the problem into a system of three tridiagonal 
matrix equations, along lines parallel to the x, y and z co-ordinate directions, which 
can be solved using the Thomas Algorithm (Godunov, 1959). Solving the three 
tridiagonal equations yields a solution for the concentration at the next time step, 
q m +i via the intermediate concentrations (both dummy variables), q* and q**. The 
convective operators are discretized using an upwind first order scheme, while the 
diffusive terms are discretized using a second order center difference scheme. The 
time operator is a simple forward difference term. The error is of 0( Ax, Ay, Az, At). 

The like terms in (3.9) are collected to obtain the following equations for a 
single spatial mesh point (n,p,r): 
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n Z » ,* 

T qm_ a * q 


(3.10) 
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(3.11) 

This system of equations is then represented as a single matrix equation in 
terms of the State transition matrices, Ai and A2. 


where 
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and 


a 2 = a' 2 j = 


0 0 ( Aj/2 + Ay + A z — r/ A t) 
0 0 -Ay/2 


0 0 — A z /2 


(3.15) 


A z , Ay, and A z are finite differenc 1 approximations used in the State 
Transition matrices representing each of the spatial directions, where, for instance, 
A x is the approximation of 


d — — d du 

di DM dx ~ di 


(3.16) 


The solution to this set of equations for the appropriate initial and boundary condi- 
tions discussed in the next section yields the concentration profiles for the cabins. 
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3.2.2 Boundary conditions 

Boundary conditions may be of the following kinds: 

Boundary Condition of the first kind: 

This is also known as a Dirichlet condition, and specifies a given concen- 
tration along a boundary. This occurs commonly at inlet ducts where the flow is 
coming in at a certain concentration of contaminant. 

Boundary condition of the second kind: 

This is also known as a Neumann condition , and specifies a concentration 
derivative normal to the surface of the boundary A wall, for example is represented 
as a Neumann condition with the derivative of the concentration being set to zero. 

Boundary condition of the third kind: 

Also known as a Robbins condition, this specifies a combination of a concen- 
tration and a flux at the boundary, and does not usually occur in cases of contaminant 
dispersion, although it occurs commonly in convection diffusion when applied to heat 
transfer problems. 

Continuity boundary condition: 

This is prescribed typically along interfaces, open boundaries, and at ducts 
linking cabins, and for no barriers to mass-transfer, specifies that the flux must be 
constant across a boundary. If a barrier exists, say, a membrane across which the 
cabin air diffuses, a resistance to mass transfer may be used to specify the boundary 
condition. 

The different boundary conditions are used in this work in the following 
cases: (Skliar, 1996; Roache, 1972) 

The nature of the boundary must be described mathematically to com- 
pletely specify the problem. This involves describing the volume of interest, as a 
single volume or as a set of volumes glued to one another, and the appropriate bound- 
ary conditions. The boundary condition can be of three main types, the Dirichlet 
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boundary condition, where the concentration at a given boundary is specified, the 
Neumann boundary condition, where the mass flux at the boundary is fixed, and the 
mixed or Robbins condition, where a function of the flux (the derivative) and the 
concentration is specified. In physical terms, the boundaries are encountered in the 
following cases. 

The wall: This is a no flux Neumann boundary condition. This also is the 
default boundary condition for the model. 

Duct In/Out: This would happen in the case of a window of some kind, 
with significant convective flow in/out. The velocity at that mesh point is then 
specified by the free stream velocity, and the diffusion is just as before. These are 
inflow and outflow boundaries and the concentration is ether allowed to float (for 
outflow ducts) or specified as a Dirichlet condition (for inflow ducts), as far the the 
mesh is concerned. The value of the float is determined by the free-stream conditions 
outside the volume and will have to be specified from time to time. 

Open Hatch: This is similar to the duct, except that flows are ignored. 

Membrane wall: This boundary is treated like an interior point, with the 
velocity set to the flow velocity determined by the membrane, and the diffusivity is 
changed to the value dictated by the membrane. 

Removal device wall: This is a specie! case, and allows for a boundary to 
exist within the volume, with a diffusivity different from the rest of the volume, 
or can act as a sink, which then would be treated a source term with a negative 
capacity. 

The solution to any partial differential equation (PDE) can depend on the 
boundary conditions and the initial conditions applied to the PDE. It is therefore, 
not surprising that the specification of the computational boundary condition, be- 
sides affecting numerical stability, affects the accuracy of the PDE solution in a 
significant manner. The intermediate values, q M and q** are not necessarily approx- 
imations to the value at the end of the iteration. As a result, particularly for high 
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order accuracy methods, the boundary conditions for the intermediate sweeps must 
be obtained in terms of the boundary values at t — nk and t = (n + l)fc. For 
boundary conditions independent of time (eg. Dirichlet conditions), the conditions 
are straightforward (LeVeque, 1985) 


q* = q** = p m+1 (3.17) 

where g m+1 is the specified concentration at the boundary. For Neumann conditions 
the relations are: 


«’ = fy* + (1 - 7-4j)(l - 7^)(^"- +I - (V) (3.18) 

9” = + (1 -jA 2 Mjg m +' - 2 9 ”) (3.19) 

where £ m+1 is the specified flux boundary condition, and 7 = (3 = Ax/2. 

The Neumann boundary also requires the calculation of the flux at the 
boundary. At the boundary, however, a center difference cannot be used to compute 
the first derivative since there are no neighbor mesh points in one of the directions 
( say, to the left of a point on the x boundary). This problem is circumvented 
by using a reflective boundary condition (where the domain is extended left of the 
actual boundary). One-sided differences could also be used. In this work, we use the 
reflective boundary condition to model the boundaries. 

3.3 Model Testing 

We first ran simulations of the contaminant dispersal using cases where we 
expected a known pattern of dispersal to see if the model was performing satisfac- 
torily. The model was therefore tested for stagnant cabins, and for cabins with air 
flow, with two kinds of contaminant sources, puffs, which are instantaneous releases 
of contaminant, and continuous streams of contaminant flow. 
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Case 1: Stagnant Room - Source in t he center. 

In the first test case, the velocity vector was assumed to be identically 
zero at all the node points. This is a useful test case, since we know a priori what 
the model should predict: a gradual spread of the contaminant outwards, while 
preserving the symmetry of the distribution. A point source was introduced at the 
geometric center of the cabin, and the contaminant concentrations were observed for 
subsequent time steps. A31x31x31 mesh was used to discretize a model cabin 
of size 7 x 7 x 7m, and Figure 3.2 shows the contours at grid point slices 1, 10, 15, 
and 29 (in the z-direction) for time Steps 3,7,29, and 49. Note that the solution is 
symmetric in all directions. 

Case 2: Room with air flow - Source in the center. 

Having tested a rather rudimentary test case, we now considered a case with 
a fictitious wind field. A wind field with u = 0.5m/s, v = 0, w = 0 was used at all the 
node points. This wind field again has the advantage that we can expect a certain 
pattern in the solution. Again, we use a 31 x 31 x 31 mesh, with two ducts on either 
side. Figure 3.3 shows the contours obtained for different time steps. Transport in 
the x-direction is mainly via convection, and diffusion in the y and z directions is 
only molecular. One can observe that the comective transport is much faster than 
the molecular transport. Consequently, we may conclude that the accuracy of the 
flow field will largely control the accuracy of the transport model. 

3.4 Contamination scenarios 

We now proceed to test the working cf our model by simulating some test 
cases. In this section, we use the wind field that we obtained in the sample cases 
(Narayan and Ramirez, 1998b) to observe hov contaminant dispersal occurs. We 
consider two separate cases, both of which represent actual Space Station contami- 
nant scenarios. 

Case 1: Steady-State Contaminat ion 
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Figure 3.2: Contaminant concentration contours for a stagnant room with a contin- 
uous release of contaminant in the center of the cabin. The contours are shown for 
4 horizontal slices per time step (left-right, top-down), at levels 0.07,0.7,1, and 1.9m 
from the floor of the cabin at time-steps 3,7,19 and 49. 
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Figure 3.3: Contaminant concentration Contours for a room with flow, with a uni- 
form velocity of 0.5 ms -1 shown at levels 0.07,0.7,1, and 1.9m from the floor of the 
cabin at time-steps 3,11,45 and 99. 
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This case attempts to mimic the routine operation of the Space Station 
module. For this specific case, we are assuming a release of 20 mg of carbon dioxide 
over the first two time steps(of 4.5s each). In addition, we have a steady input 
of carbon dioxide in the inlet streams. The cabin has a residual carbon dioxide 
concentration of 0.3 volume % for its initial condition. The THC air comes at a 
concentration of 0.5 volume % and the IMV with a concentration of 0.71 volume %. 
This would be realistic since the THC is treated , and we could assume that one of 
the other modules has a great deal of astronaut activity, and thus the high carbon 
dioxide level. The simulation was done for about 300 time steps (about 1350 s). By 
this time, we observe a steady-state concentration distribution. Figure 3.4 shows the 
contamination levels at four different slices in the cabin, at levels 0.067m, 0.67m, lm, 
and 1.87m from the floor of the cabin after about 1300 s. The origin, (0,0,0) refers 
to the point in the cabin at the left bottom corner in Fig. 2.1. The surface plots of 
the slices closer to the middle of the room show similar profiles, and the exits and 
the consequent drop in concentration levels of the contaminant are clearly visible. 
This is due to the strong blast (relative to the rest of the cabin) of wind removing 
the contaminant through convective transport. 

Case 2: Sudden Release of Carbon dioxide 

Carbon dioxide may be used to extinguish a fire. A large release of carbon 
dioxide then would occur over a small time frame, and we wish to monitor how the 
contaminant levels gradually decrease. Figures 3.5 and 3.6 show the concentration 
levels 90 and 1500 s after the release. The surface plots and contours are shown 
at planar slices 0.067m, 0.67m, lm, and 1.87m from the floor of the cabin. During 
these 1500 s, more than 70% of the released carbon dioxide has been flushed out from 
the room. Figure 3.5 shows the profiles 90 s after the release, which happened near 
the bottom left corner of the cabin (0.27m,0. 14m, 0.17m). In this figure, substantial 
amounts of the carbon dioxide are still present near the location of the occurrence, 
although the levels drop off near the outlets. In Fig. 3.6, which shows the profiles 




Figure 3.4: Steady State Contamination CO 2 profiles for the cabin with an initial 
carbon dioxide concentration of 0.3 volume %, an IMV inlet flow at 0.71 volume %, 
and a THC concentration of 0.5 volume %. The concentration profiles are shown at 
four levels, 0.067m (level = 1), 0.67m (level = 1)), lm (level = 15), and 1.87m (level 
= 28) above the floor of the cabin. 
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1500 s after the release, we observe that the carbon dioxide is almost well-mixed with 
only a slight bulge at the location of the release in the top slice. This accumulation 
of carbon dioxide near the ceiling is a consequence of the low flow region at the top, 
and because of the stagnation zone that exists near the ceiling. By this time, it is 
the flow profile which largely decides the contamination concentration profile, which 
is indicative of the fact that convective transport is the dominant mass transfer 
mechanism. For a sample average flow velocity of 1.5 m/s, the mass transfer Peclet 
number is 2 x 10 5 , which implies that convective transport dominates molecular 
diffusion. 

Note here that there is a significant variation in the concentrations across 
the room, which might mean that lumped models of the cabin would be grossly 
inaccurate. Also, there are regions of accumulation in the room. This could mean 
that SMACs could be locally violated, even though the concentration averaged over 
the entire cabin may be below the SMAC limit. The flow field is an important 
parameter in the way contaminants spread through the cabin and needs to be closely 
monitored. 

3.5 Modeling of cleanup process 

Crucial to the Advanced Environmental Monitoring process is the modeling 
and monitoring of the contaminant removal processes. Removal processes range from 
HEPA filters, molecular sieves, and adsorption packs to the use of the Sabatier and 
the Bosch processes for the regeneration of carbon dioxide. 

The approach used in this work is that all removal devices can be modeled 
in two distinct ways; either as devices that remove a certain percentage of the con- 
taminant and leave the rest in the medium, or as devices that remove almost all of 
the contaminant leaving only a small residual concentration behind. The air that 
passes through the removal device is then assumed to return to the main stream and 
re-enter the module. An additional parameter is the amount of contaminant that 





Figure 3.5: CO 2 profiles at four different elevations, about 90 seconds after a sud- 
den release of carbon dioxide at grid location (4,5,8) with an initial carbon dioxide 
concentration of 0.3 volume %, an IMV inlet flow at 0.71 volume %, and a THC 
concentration of 0.5 volume %. The concentration profiles are shown at four levels, 
0.067m (level = 1), 0.67m (level = 10), lm (level = 15), and 1.87m (level = 28) 
above the floor of the cabin. 


Figure 3.6: CO 2 profiles at four different elevations, about 1500 s after the sudden 
release of carbon dioxide at (4,5,8) with an initial carbon dioxide concentration of 
0.3 volume %, an IMV inlet flow at 0.71 volume %, and a THC concentration of 0.5 
volume %. The concentration profiles are shown at four levels, 0.067m (level = 1), 
0.67m (level = 10), 1m (level = 15), and 1.87m (level = 28) above the floor of the 

C£ 
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Table 3.1: Samp! 


ie parameters used for the removal devices 


Capacity 

Boundary-type 

Removal Efficiency 

Name 

100 

1 (percentage ) 

0.95 

C02 scrubber 

50 

2 (Dirichlet) 

0.0004 

adsorber 


the removal device can safely accumulate before it stops functioning. 

In this research, the removal of all substances is monitored by a subroutine 
that calculates the flux of each substance at all the outlets, monitors their removal 
based on the input parameters, and updates the new inlet concentration boundary 
conditions. A sample simulation, with the sample parameters of Table 3.1 was carried 
out using a cabin with the removal devices in place. A source of 60 mg/L was emitted 
over a 10 second period, and the THC and IMV inlets had concentrations of 0.0015 
mg/L and 0.002 mg/L, respectively. The concentration profiles after 450 s are shown 
in Figure 3.7. 


This subroutine is flexible enough to model a variety of devices, including the Carbon 
Dioxide Removal Assembly (CDRA) that is part of the current baseline technologies 
in the Air Revitalization. The CDRA prototype tested(Barker et al., 1991) used a 
carbon dioxide removal rate that was a function of the inlet carbon dioxide partial 
pressure, PPco 2 an ^ was represented by the following equation: 

Removalrate(lb/hr) = 0.1579 * PPco 2 ( mm #<?) ~ 0.0348 (3.20) 

The equation is valid for carbon dioxide partial pressures between 2.0 and 3.9 mm 
Hg, and can be easily included in our removal \ ubroutine. 

3.6 Summary and Conclusions 

In this Chapter, the transport model was developed and implemented for 
air contaminant dispersal aboard spacecraft. The transport model uses the flow field 




Figure 3.7: Simulation of the cabin concentrations with the removal devices in use. 
The parameters of the removal devices is shown in Table 3.1 and are placed at the 
IMV and THC outlets. The concentration surfaces are shown at four horizontal 
slices in the cabin, at vertical levels 1,10,15, and 29 (grid location in the z-direction 
out of a total of 30). 
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calculated using methods described in Chapter 2 as an input along with the cabin 
geometry and the fluid parameters, and predicts the dispersal of air contaminants as a 
function of time. This transport model is a three-dimensional model, an improvement 
over previously developed lumped and two-dimensional models, which is shown to 
be useful when there are wide fluctuations in the contaminant concentration in the 
cabin such that local violations of the SMAC level can occur even when the average 
across the cabin stays below the SMAC level. Convective transport was found to 
dominate molecular diffusion, evidenced by a mass transfer Peclet number around 2 
x 10 5 . The transport model accounts for sources of the contaminant, both in the inlet 
flows and from inside the cabin, and is also extended to model the effect of removal 
devices that are commonly used in the Space Station. The model is accurate for 
monitoring purposes, and is also computationa ly suitable for real-time application. 



Chapter 4 


State estimation using Implicit Kalman filtering 


4.1 Why estimate the state? 

State estimation is necessary while monitoring any physical process because 
there are always uncertainties-faults in the process under observation, errors in the 
mathematical model that is assumed to adequately describe the given physical pro- 
cess, and changes in parameters that can cause real concentrations to be different 
from those predicted by the model. The objective of the filtering process is to arrive 
at an estimate that is unbiased, i.e. has the smallest error in the least-square sense, 
and which gives one an accurate picture of the actual system. The cost of sensors is 
high, both in terms of the monetary expense, and on account of weight and electrical 
power issues, which necessarily restricts the number of sensors that can be carried 
aboard. This gives rise to the issues of placement and selection of sensors, which is 
an area of active research (Smith, 1996). The estimation process is very crucial to 
the fault detection and diagnosis process since the matrices and calculations used 
in fault and diagnosis procedures are used to make inferences about if and where a 
fault (contaminant leak) has occurred. 

An identification process usually has at least three main ingredients: 

• A priori knowledge in the form of a mathematical model about the unknown 
system and the noise. 

• A measurement system that provides discrete or continuous measurements 
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of a subset of the state. 

• A performance measure of the identification algorithm. 

The state estimation procedure closely follows that proposed by Skliar and 
Ramirez (Skliar and Ramirez, 1997b). The estimation problem, formulated simply 
is as follows: Given a stochastic process that represents a dynamic system, we are 
interested in knowing the value of x(k) for some fixed fc, where x(k) is not directly 
accessible to us for observation. We have a sequence of measurements that are 
causally related to x(k) by means of a measurement system M and measurement 
data z(z), and we wish to utilize these data to infer the value of x{k). We denote 
the estimate of x(k) by x(k) and define it to be some n-dimensional, vector-valued 
function 0* of the measurements, viz., 

x{k\j) = <f>k[z(i),i = (4.1) 

where k refers to the time when the estimation is made, j refers to the time until which 
measurements are taken and used, and i is an index that refers to the measurement 
signal being used. Crucial to the estimation p ocess is the definition and notion of 
the estimation error which is defined by the re at ion 

£(%') = x(k) - x(k\j) (4.2) 

Ideally, x = 0 and the estimate is exact. When this is not the case, we 
assign a penalty for the incorrect estimate. Th s is done through a penalty or a loss 
function L which has the following properties: (Meditch, 1969) 

1. The loss function is a scalar-valued function of n variables. 

2. L(x = 0) = 0. There is no penalty if the estimate is exact. 

3. L is a non-decreasing function of the distance of the error from the origin 


in n-dimensional Euclidean space. 
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4. The loss function is symmetric about the origin in the same n-dimensional 
Euclidean space. 

4.2 Implicit Kalman Filter 

One of the classical methods of state estimation is the well established 
Kalman filtering algorithm (Kalman, 1960; Kalman and Bucy, 1961). The Kalman 
filter, which has many different implementations now, and which is widely used for 
the purpose of state estimation for dynamic systems that have random perturbations, 
is an unbiased and minimum error variance recursive algorithm to optimally estimate 
the unknown state of a dynamic system from noisy data taken at discrete real- 
time. In the Kalman filtering paradigm, the uncertainties of the model and the 
measurements are represented by additive stochastic white noise. 

In addition, the measurements, z(i ) and the measurement errors e(z) are 
assumed to possess the following properties: 

(1) The measurement errors have a zero mean, 

E(e(i)) = 0, 

where E is the “expected value operator”. A zero mean error is different 
from a random error in that a random error may have a non-zero expected 
value, which is also known as a bias. Here, we are assuming that our sensors 
are not biased in any particular direction. 

(2) The measurement error has a constant variance, which does not change with 
time or with any other parameter 

(3) The errors are additive, i.e, 


z(i) = x(i) 4- e(i) 
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(4) The measurement errors are uncorrelated, i.e. for two measurement errors, 
c(o) and e(p), 

cov(e{o),e(p)) = E{[e(o) - £(e(o))][c (p) - £(e(p))]} = 0 for o ^ p (4.3) 


The Implicit Kalman filter is one of the many alternative Kalman filters that 
have been developed over the years, which are theoretically equivalent to the original 
formulation, and which have been formulated for various desirable features (Carl- 
son, 1990; Chin et al., 1995; Jordan, 1967) including enhanced numerical stability, 
computational accuracy, reduced computational requirements or for implementation 
using parallel computing (Jover and Kailath, 1986; Morf and Kailath, 1975; Paige 
and Saunders, 1977; Roy et al., 1991). Many of these variations are also discussed 
in textbooks (Chui and Chen, 1991) dealing with the topic of Kalman filtering. The 
Implicit Kalman Filter was shown to be particularly efficient in treating descriptor 
systems with sparse transformation matrices, an example of which is the convection- 
diffusion equation that forms the core of our mathematical model for air contaminant 
dispersion. 

If we recast the model equations for the transport model developed in Chap- 
ter 3 and include the additive noise, the model can be written as a single matrix 
equation 


where 



fm 


C(m) 

AlQm+l = A2Q m + 

0 

+ 

0 
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Wim 


Qm 
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** 

m 


Qm 


(4.4) 


(4.5) 
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C(m) represents the stochastic disturbance transition matrix that acts upon wi, 
which represents a stochastic disturbance that is an uncorrelated Gaussian white 
sequence with a zero mean and 


£[wi m w lm T ] = Q 


(4.6) 


Q is a diagonal matrix that represents the model noise, and is closely tied to the 
uncertainty surrounding the mathematical model for our physical process. A highly 
accurate model would have low values for its Q, and one that does not represent the 
physical process too accurately would have high values for its Q. 


and 


and 


A, = {Aj J } 


( — A r 0 

v 2 At) u 

JL ( ~ A y E_ \ 

At v 2 At) 
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r 

M 



A 2 = A-2 = 


0 0 (A^/2 + A y + A z — v j A t) 
0 0 -A y /2 

0 0 —A z /2 


z m+ 1 


0 0 H(m + 1) 


Qm+l + v m+l 


(4.7) 


(4.8) 


(4.9) 


v m+ i is an uncorrelated Gaussian white noise sequence that represents the mea- 
surement noise with a covariance represented by R. The R matrix is a diagonal 
matrix that contains information about the uncertainty surrounding the measure- 
ments. Lower values for R indicate more accurate sensors. The estimation of the 
contaminant concentration is determined from the sequential solution to the follow- 
ing tridiagonal equations: 
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A x r 

2 At q 

— ( 2 ^ y ^ 2 + f m 

(4.10) 


+Li[z - Hiy m _f_ X | m ] 


, A y r 

A v r * _ . ^ . 




(4.11) 

(_ T - At )q ™ +1 

A. 2 r 

= — ^qm - -^q** + l 3 [z - Hiy m+1 | m ] 

(4.12) 


where the predicted estimation of the auxiliary variable y — A2Q is given by the 
following equation: 

Ai 3 [ f 

a 2 1 m 

ym+l|m A^ 3 4m|m 0 

A^ 3 0 

The superscripts in the matrices A 2 3 , A 2 3 , and A 33 refer to the row and column 
partitions of the A2 State transition matrix. 

Note that the Eq. 4.10-4.12 axe identical in structure to the model equations 
in 3.11, the only difference being the addition of the last term, which is the effect of 
the filter on the equations. 

The modified measurement matrix Hi = {Hij} is calculated using the 
following equation: 

Hu H 12 H 13 Ai= 0 0 H , (4.14) 

which is equivalent to solving the following equations: 

H 13 A 33 = H, 


(4.13) 


H 12 Af 


-h 13 a 32 , 


(4.15) 



65 


HnA} 1 = -H 12 Af. 


The Implicit Kalman Gain, a matrix that multiples the estimated error residual is 
given by the equation 


Wm = PIUm.XIH.P' + 1|m Hf + R (m + 1)]-' (4.16) 


where Pf n+1 | m is the predicted error covariance matrix given by the equa- 


tion 


P y 

m+l|m 


A 2 3p llm A 2 3T 

^ m\m * 

+cqc t 

a^p^a^ 

AfP^|„Aj 3T 


A| 3 PL |m Af r 

A?P^ |m Af r 

A?Pl, m Af T 


A 2 3 Pl| m A 2 3T 

z m\m * 

AfP 9 Af r 

* m\m * 

AfP’ , A^ 

* m\m £ 


(4.17) 


The actual estimation error covariance matrix, P y m +i| m +i is then calcu- 
lated from the predicted error covariance matrix using the relation 


P y m+l|m+l — [I “ Pm+lHi]P y m+1 | m , (4.18) 

The error covariance matrices are measures of the uncertainty inherent in the com- 
puted quantities. 

The error covariance in the state, P^ +1 | m+1 which is the uncertainty asso- 
ciated with the computed concentration is determined from the equation 

P y m+l|m+l = A!P^ +1 | m+1 Af (4.19) 

The solution to the last equation is reduced to a sequential solution of the 
following six tridiagonal equations: 
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where [ P m+l|m+l] 33 = P m+l|m+r 

We can now formulate the algorithm of the estimation of the contaminant 
concentration q m +i based on the measurement data and the transport model. 

1. Compute the predicted estimate of the concentration, y m +i| m by prop- 
agating the concentration at the previous time step, q m | m according to Eq. 4.13. 

2. Successively solve three tridiagonal matrix equations for the modified 
measurement matrix U\. 

3. Solve the tridiagonal model equations with the new perturbation to 
obtain the optimal estimate q m +i|m-f-i* 

The calculation of the gain L m+1 of the Implicit Kalman Filter follows the 
following algorithm: 

1. Calculate P^ +1 | m according to Eq« 4.17. 

2. Calculate the Implicit Kalman filter gain from Eq. 4.16. 

3. Calculate P^n +1 | m+1 according to Eq. 4.18. 
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4. The final step is to initiate the gain calculation for the next time step, 
for which one sequentially solves Eq. 4.21. This is by far the most time consuming 
step. The updating step can also be implemented through the use of square-root 
filtering schemes, which can reduce the computation time, and also provide some 
added stability to the filter (Skliar and Ramirez, 1997b). 

4*3 Filter Implementation and Testing 

A “true” test of the filter can only occur in an experimental setting, with a 
physical cabin and measurements. In the absence of that, we tested the filter using 
the results of the model itself. We added a random Gaussian noise to the contaminant 
concentrations from the model, as we would expect to get in a real setting and then 
checked to see if the filter was able to track the contaminant concentrations with 
sufficient accuracy. 

The filter is a computationally intensive program, and we found that it was 
taking around 30 seconds for each time step on a DEC-Alpha Station 250 4/266. 
We therefore ran the model and the filter with a 30s time step. We are using a 
set of five sensors to estimate the concentration in the whole room. In order to 
test the filtering process, we have crowded all of the sensors in one corner of the 
cabin, to see how robust the filtering process is. Table 4.1 shows the position of 
the sensors. The co-ordinates refer to the location of the sensor grid points within 
a rectangular geometry of size 6 x 2 x 2m which used a grid of dimensions 15 x 
28 x 30. For the filter to be useful, good estimates of the model uncertainty and 
the measurement uncertainty are needed. This estimation is not a trivial process, 
and will require among other things, experience in running the filter for specific 
systems and conditions. The model uncertainty, Qdiag can be closely tied to the 
amount of turbulence in the system and the error in the numerical solution to the 
convection-diffusion equation. The measurement uncertainty, Udiag should reflect 
the uncertainty present in each sensor, which is quantifiable by test experiments. 
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Table 4.1: Sensor location and the associated measurement noise 


Sensor 

Co-ordinates 

R 

i 

3,2,3 

0.001 

2 

1,1,3 

0.001 

3 

2,2,2 

0.001 

4 

3,1,2 

0.0019 

5 

2,4,3 

0.0019 


Estimation of Contaminants 



Figure 4.1: Filter Performance-Tracking an arbitrary point in the cabin. The real 
concentration at point (1,2,3) is given by the solid line, and the estimated concen- 
tration from the filter is given by the dotted line. 
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Figure 4.3: Filter Performance-TYacking at Sersor #1, located at (3,2,3). The solid, 
dotted, and dashed lines indicate the true concentration, the measured concentra- 
tions, and the optimal estimate of the concentrations using the Implicit Kalman 
Filter, respectively. 
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0 100 200 300 400 500 

time 


Figure 4.4: Filter Performance-Tracking at Sensor #2, located at (1,1,3). The solid, 
dotted, and dashed lines indicate the true concentration, the measured concentra- 
tions, and the optimal estimate of the concentrations using the Implicit Kalman 
Filter, respectively. 




72 


The graphs show sample results of the contaminant concentration estima- 
tion for 30 time steps. Figures 4.1 and 4.2 show the tracking of two arbitrary points 
in the domain, while Fig. 4. 3-4.4 show how the filter functions at the sensor loca- 
tions in the presence of measurement noise. The tracking is fairly accurate and is 
mathematically consistent. The error for points in the domain that are distant from 
the sensors is naturally higher. The filter slightly lags the “true” values. That, too, 
should be expected, since the filter does not “know” about new source emissions. 

Figure 4.5 shows the error at four different sensor locations as a function 
of time. The dashed line, which is the bound of the error is determined by the 
estimation error covariance at that location. In this example, we have chosen 3 a 
as the bound, where a represents the standard deviation of the expected estimation 
error, based on the fact that 99 % of the estimates would fall within this bound. It is 
expected that the error stays within the bounds that are indicated by the uncertainty, 
which changes from its initial value until it reaches a steady state value. Figure 4.6 
shows how the uncertainty varies with time. 

4.4 Uncertainty 

4.4.1 Modeling uncertainty 

The main source of uncertainty in the system being studied is that of the 
wind velocities. The flow field is by far the mos , important parameter in the disper- 
sal of the contaminants. In addition, there axe t le other usual uncertainties inherent 
in the physical modeling of any system, errors in the measurement, and the pres- 
ence of faults. There are two reasons why there are uncertainties in the flow field. 
Firstly, turbulence is stochastic in nature, and the flow field obtained, by definition, 
is approximate since it is a time averaged quant ty. Graphs showing the actual mea- 
sured velocities indicate this quite clearly (Zhang et al., 1992; White, 1974). The 
other reason is that numerical procedures essentially yield approximate results. The 
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Figure 4.5: Filtering results showing the noisy measurements and the filtered es- 
timates; the ordinate represents the measurement noise (dotted line) and the es- 
timation error (solid line) as a function of time. The dashed line represents the 
error bounds, and is equal to 3cr, where a represents the standard deviation of the 
expected estimation error 
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residual norms of the errors are a measure of how approximate the numbers are. 
The objective in studying the uncertainty is to structure the uncertainty in order to 
produce a measure of how reliable our estimate of the state is and how robust the 
procedure is. It is crucial that the diagnosis be accurate, that no fault go undetected; 
at the same time, any false alarms must be avoided since that would mean the loss 
of precious crew time. Most fault detection procedures will ultimately detect any 
given fault; the objective is to detect the fault as early as possible using a proce- 
dure that is known not to yield too many (or any) false alarms. The structuring of 
the uncertainty therefore, is a prerequisite to the fine tuning of the fault detection 
algorithm. 

An analysis of the uncertainty would require an experimental set-up, and 
data in order to evaluate the model results. In the absence of that, one could 
artificially change some parameters, and then evaluate the performance of the filter 
and the model. We next consider the effect of a change in the inlet velocity on the 
model results. 

4.4.2 Effect of randomness in inlet velocity 

The dispersion model assumes the existence of a steady state flow profile. 
Here, we examine the effects of a variation in the inlet velocity on the estimation 
error. 

For the turbulent flow field simulated in Chapter 2, the THC inlet velocity 
is specified to be 0.5 m/s. Now, assume that the velocity decreases by 15 %. 

The steady-state Navier-Stokes equations are next solved for the new bound- 
ary conditions. A study of the actual velocities shows that some internal velocities 
decrease by more than 15 %, depending on the location, not a surprising result since 
the Navier-Stokes equations are non-linear, and the turbulent energy itself causes 
substantial noise in the system. 

The new velocities are used in conjunction with the mathematical model 
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Table 4.2: Sens or location and the associated me asurement, noise 

t — 1 1 1 


Sensor 

Coordinates 

R 

0 

10 10 10 

0.0001 

1 

8 2 23 

0.001 

2 

3 5 7 

0.0001 

3 

6 6 6 

0.0023 

4 

4 11 22 

0.019 

5 

8 14 18 

0.0019 


to generate the “true” concentrations for the given concentration initial and bound- 
ary conditions. The filter, however, is not updated with these velocities, since our 
objective is to evaluate its performance when conditions change unknown to the 
model. Monitoring the inlet flow velocity would take care this change and update 
the model, a procedure that is discussed in the next section. In the first case, we use 
the uncertainties from Table 4.2. 

The tracking at Sensor locations 0 and 4 is shown in Figs. 4.7 and 4.8. 
The filter uncertainties in this case balance the model uncertainties, and the filter 
consistently over-predicts the concentrations at both locations. The main reason 
for this is the fact that the reduced inlet velocities causes a reduced flux of the 
contaminant into the chamber, since the inlet concentrations remains unchanged. 
This 15 % reduction in the flux persists throughout the duration of the experiment, 
and the filter, constrained by the mass balance over predicts the concentration. 

The residual error curves for Sensors 0 and 1, shown in Figs. 4.9 and 4.10 
show the error that is negative (with the estimated concentrations being lower than 
the actual) and increases with time. The errors increase at different rates, and the 
rates depend on the local velocity and location with respect to the ducts. 

The Euclidean norm of the prediction error, shown in Fig. 4.11 rises dra- 
matically, and clearly exceeds the error bounds, and the fault in the system is quite 
apparent. 

This case dealt with the situation in w lich the model and sensor uncertain- 
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Figure 4.7: Tracking of the contaminant concentration at sensor location #0 with a 
15 % reduction in the inlet velocity. The filtered estimate (dashed line) over-predicts 
the real concentration (solid line). The measurements are shown by the dotted line. 
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Figure 4.8: Tracking at location 4 with a 15 % reduction in the inlet velocity. The 
filtered estimate (dashed line) over-predicts th< real concentration (solid line). The 
measurements are shown by the dotted line. 
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Estimation error for sensor- 0 



Figure 4.9: Residual error at sensor location 0 with a 15 % reduction in the inlet 
velocity. The magnitude of the error is increasing with time. 


Estimation error for sensor— 1 



Figure 4.10: Residual error at sensor location 1 with a 15 % reduction in the inlet 
velocity. The magnitude of the residual (estimation) error increases with time. 
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Euclidean norm of the estimation error 



Figure 4.11: The Euclidean norm of the erroi with a 15 % reduction in the inlet 
velocity is shown increasing with time. The estimation error clearly exceeds the 
error bounds (dashed line), indicating the presence of a fault. 
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Table 4.3: Sensor location and the associated measurement noise- very low sensor 
uncertainty 


Sensor 

Co-ordinates 

R 

0 

10 10 10 

0.000001 

1 

8 2 23 

0.00001 

2 

35 7 

0.000001 

3 

6 66 

0.000023 

4 

4 11 22 

0.00019 

5 

8 14 18 

0.000019 


ties were weighted equally. If we expect extensive noise in the system, for example, 
if a 15 % change in the inlet velocity were common, one would then have to give the 
measurements a higher certainty level in order to obtain better tracking performance. 

For this situation, we repeat the previous case, but with the lower mea- 
surement uncertainties listed in Table 4.3. The sensor uncertainties, Rdiag have been 
reduced to about one hundredth of their usual values. The results are shown graph- 
ically in Figs. 4.12 through 4.16. One can clearly observe that the tracking using 
these very low uncertainties is an improved over the previous case. In this scenario, 
the model is relatively unimportant, and the measurements become paramount. This 
all comes at a cost, of course. With these very low sensor uncertainties the ability 
to filter out the noise is impaired. 

4.4.3 Double Filter 

Since one of the major sources of uncertainty and error in this model is the 
velocity, a variant of the Implicit Kalman filter such as a double filter could be used 
for improved performance. The double filter is a Kalman filter, but the difference 
between the measured and model velocity is also used to update the final state of 
the system. It is expected that this will increase the sensitivity of the filter and 
consequently increase accuracy. The filter will therefore have a set of concentration 
sensors, measurements from which are used to update the right side of the state 
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Figure 4.12: Tracking at sensor location 0 with a 15 % reduction in the inlet velocity 
for very low measurement uncertainties. The filter (dashed line) now tracks the 
measurements (dotted line) very closely. 
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Figure 4.13: Tracking at sensor location 4 with a 15 % reduction in the inlet velocity 
for very low measurement uncertainties. The filter (dashed line) still over-predicts 
the concentration, but the tracking is slightly better than in Fig. 4.8. 
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Estimation error for sensor— 0 



Figure 4.14: Residual error at sensor location 0 with a 15 % reduction in the inlet 
velocity for very low measurement uncertainties. The presence of the residual error 
indicates the presence of a fault, but the error in this case is lower in magnitude than 
with higher sensor uncertainties. 
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Estimation error for sensor- 1 



Figure 4.15: Residual error at sensor location 1 with a 15 % reduction in the inlet 
velocity for very low measurement uncertainties. The presence of the residual error 
indicates the presence of a fault, but the error in this case is lower in magnitude than 
with higher sensor uncertainties. 
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Euclidean norm of the estimation error 



time (s) 


Figure 4.16: Euclidean norm of the error with a 15 % reduction in the inlet velocity 
for very low measurement uncertainties. The presence of a fault in the system is 
evidenced by the fact that the residual error solid line) clearly exceeds the error 
bounds (dashed line). 
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estimation equations, and a set of velocity sensors, which will perturb the left-hand 
side of the state estimation equations. The new state estimation equations, will 
therefore be of the form, 


(-a«/2 - 1/A t)q* 
(-«»/ 2 - 1/A t)q** 
~ l/At)q m+ i 


( a x /2 + cty + ot z 1/A t)q m + / m + Li[z — 77iy m -f-i| m ] 
—a y /2q m — l/Atq + £ 2(2 — #iy m +i| m ], (4-21) 

—^-Qm ~ l/Atq** + Lz[z - //iy m +i| m ], 


where a x , a y , and a z represent the filtered left hand side terms. The time step used 
in updating the velocity filter could be much larger than the time step used in the 
rest of the filter, in order to minimize the computations needed. 

Implementation of this double filter was out of scope for this work, but can 
be easily implemented and tested. 

4.5 Summary and Conclusions 

A state estimation procedure was implemented using the Implicit Kalman 
Filter, which provides accurate estimates of the contaminant concentrations at all 
points in the cabin, using the transport model developed in Chapter 3 in conjunc- 
tion with a measurement system. The filter is an effective tool for rejecting sensor 
noise, and provides smooth estimates of the state of the system in real time. The 
performance of the filter in the presence of a major disturbance was studied, which 
showed that a proper choice of the model and measurement noise covariance matrices 
can lead to good tracking behavior in the presence of noise. Parameter estimation, 
and its use in combination with the state estimation procedure could lead to further 
improvements in overall estimation performance. 



Chapter 5 


Fault detection in distributed parameter systems 


5.1 Introduction 

In this chapter, we discuss the use of the Implicit Kalman Filter in the 
implementation of a Fault Detection algorithm. Fault detection is the procedure 
which alerts the user to a malfunction in the system. Fault detection is the first step 
in the comprehensive Fault Detection and Isolation (FDI) problem. Basseville, in his 
discussion of current methods in FDI (Basseville, 1997), mentions that FDI is split 
into two steps; the generation of residuals , which are ideally zero under fault-free 
conditions, minimally sensitive to noises and disturbances, while being maximally 
sensitive to faults, and residual evaluation , wiich concerns the design of decision 
rules based on these residuals. A detailed account of fault detection algorithms can 
be found elsewhere (Basseville and Nikiforov, 1993). Methods of fault detection are 
classified into methods that axe Model based or those that are Statistically based. 
Statistical methods are the only options whe l detailed model information is not 
available and consists of continuously examining the statistical properties of the 
measurement data, and noting any substantial deviations from a pre-determined 
threshold band. Model based methods, on the other hand, use the knowledge about 
the system and infer unmeasurable characteristics of the system from the measurable 
using this knowledge. 

The fault detection algorithm consists of a sensor fault test procedure and 
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a process fault detection procedure. The sensor fault testing procedure forms the 
inner shell of the algorithm, since the sensor readings will have to be validated before 
further processing. Then, the readings are evaluated for a possible process fault. 
Under normal conditions, both the tests would be negative, and the filter would 
continue onto its next time step. For the purposes of air contaminant monitoring in 
spacecraft, we envision two kinds of faults; Instrument or sensor faults and Process 
faults . A Sensor fault, as its name suggests, implies that one or more of the sensors 
is not functioning. The fault could be in the form of a total malfunction, where the 
sensor readings are totally random with no physical basis, or could be manifest as a 
bias in some direction. 

A process fault may be present as an unknown source problem, or as a 
violation of a safety requirement. While we expect that faults will be infrequent, a 
fault detection procedure is crucial since it is under conditions of a fault that the 
utility of the system is realized. 

5.2 Sensor Fault 

One method for detecting sensor faults is that of hardware redundancy 
(Emami-Nacini et ah, 1986), in which sensors are used at each location, and an 
agreement between all three sensor readings is necessary for that sensor reading 
to be accepted as valid. If one of the three sensors shows a deviation statistically 
significant from the other two, then that sensor is considered to be faulty. Hardware 
redundancy is costly because of its need for triple the amount of hardware especially 
in a space environment, where weight are power requirements are concerns. Hardware 
redundancy, though might become feasible for on-board applications if the cost and 
weight of the sensors drop enough to offset the cost of the extra computational 
requirement due to analytical redundancy. 

Another approach is called analytical redundancy. Analytical (or func- 
tional) redundancy is a model-based fault detection procedure where a single set of 
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measurements is used in conjunction with a model and a detection algorithm in order 
to detect and distinguish between faults. Here, the set of sensors is separated into 
two. In our example, Sensors 1-3 (Table 4.1) sure part of Filter bank 1; Sensors 4-6 
part of Filter bank 2, and the complete set of sensors 1-6 form part of Filter Bank 0. 
The three filters are run simultaneously, one with the complete set of sensors, and 
one with each individual set of sensors. Three residuals characterize these banks of 
sensors, 


r o z m + 1 Hiq m 

(5.1) 

r x = z m+l - Hl^m Klm+l 

(5.2) 

r 2 = Zm+l " H l 2 qm+l|m+l 

(5.3) 


where the superscripts refer to the appropriate sensor banks. When all three banks 
show a similar residual, the sensor system is working normally. When two of the 
banks deviate, then one of the sensors in that bank is malfunctioning, and the other 
bank alone should be used in the filtering process. Mathematically, this amounts to 
checking if 

Ikolh =|| ri || 2 = | r 2 II 2 (5.4) 

within bounds of error at each time step. |jr|| 2 ^presents the Euclidean norm of the 
quantity, which represents the distance of the vector from the origin. 

Figure 5.1 shows the residuals calculated for three different banks in the 
absence of a fault, while Fig. 5.2 shows them in the presence of a fault in Bank 0 
and Bank 1. In Fig. 5.1, the residuals of the estimation error in all three banks is 
around 0.015 throughout the time period undei consideration. In Fig. 5.2, however, 
while Bank 2 maintains its residual of 0.015, Banks 0 and 1 have residuals around 
1.5, which are significantly higher. Bank 2 is therefore functioning normally, and its 
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estimates alone should be used to monitor the contaminants until the faulty sensor 
is identified and rectified. 

This method requires high computational power but does not need an in- 
terruption of the operation once a fault is detected, since a smooth transition to the 
functioning bank of sensors can be made. Alternatively, one can just monitor each 
individual sensor for its residual and determine a malfunction. Figure 5.3 shows the 
residual error for three sample sensors in the case of a sensor fault. It is quite ap- 
parent from the three sensor error graphs that Sensor 2 is malfunctioning, and that 
the others are functioning normally. The working of the filter can be seen. Since 
the filter does not “know” that a fault has occurred, the filter is starting to respond, 
and even the non-faulty sensors are showing errors that are almost exceeding the 
bounds. 

This procedure required less computational power, but the filtering process 
after the detection of the fault becomes complicated, because of the process of taking 
out the malfunctioning sensor’s readings from the existing filter. 

5.3 Unknown source 

The problem of an unknown source could be something as minor as higher 
than expected carbon dioxide levels because of increased activity in the cabin or 
could be a major leak. An unknown source will cause the model to vary consid- 
erably from the measurements and cause large estimation errors. Monitoring the 
estimation error is the key to identifying the presence of unknown source substances 
and contaminants. The way this is detected is through having safe bounds for the 
residual error, and if the residual error exceeds the bounds, then the diagnosis for 
an unknown source is initiated. Figure 5.4 shews a sample result where the error 
bound is exceeded owing to an unknown source. The case of an unknown source is 
discussed in greater detail in the next section, in a simulated Space Station scenario. 



Estimation error Estimation error 
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Euclidean norm of the estimotion error — Bonk 0 



time (s) 


Euclidean norm of the estimation error— Bank 1 



time (s) 


Euclidean norm of the estimation error - Bank 2 



Figure 5.2: Residuals of the three Filter Banks with a fault in Sensor #2. Note 
the huge change in the estimation error in the faulty banks (0 and 1) from normal 
operation (Bank 2). 
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Residuol error ot Sensor location § 2 



time (s) 



Figure 5.3: Residual error in three sample sense rs with a fault in Sensor #2. Sensor 
#2 has residuals that are clearly outside the error bounds, indicative of a sensor 
fault. The fault in Sensor #2 is beginning to affect the other sensors as well. 
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5.4 Operation fault 

An operation fault is a rather serious situation that arises when the space 
operation itself is in identifiable danger because the space atmosphere is seriously 
contaminated. One way of recognizing a fault in the system would be by studying 
the state of the system, which in this case would be the concentrations of each 
contaminant at every grid point. A fault can be posed using the norms of these data, 
the way it is commonly done in linear system theory. For purposes of contaminant 
monitoring aboard the Space Station, a fault can occur both in the I 2 and the l\ 
sense. A fault in the l\ sense occurs for the case of substances that cannot exceed a 
certain SMAC, and which is an acute toxic. For such substances, a fault occurs if 

IIM<7)II>9c ( 0.5) 

where q c is the appropriate SMAC for that contaminant, and l\ and 1 2 are the 1 and 
2 norms of the concentration vector. A fault in the I 2 sense occurs when the process 
fault affects the health of the cabin in an overall sense. For example, a fault in the 
oxygen system occurs if 

II h(q) ll< i.c 

and a fault in the carbon dioxide system occurs if 

II h{q) II > < c 

Situations could vary from high CO^ levels that have exceeded the long 
term or short term SMACS, or with toxic releases detected at levels that are known 
to be harmful to humans aboard. The situation could further be subdivided into 
either a local or a global fault, and the detection algorithm is able to distinguish 
among these. 
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5.5 Scenarios 
5.5.1 CO 2 operation 

We will discuss the utility of the monitoring system in two specific scenarios, 
both relating to the CO 2 levels in the Space Station (Narayan and Ramirez, 1998a). 
According to NASA standards, CO 2 must be present in the cabin at levels between 
0.3 - 0.80 volume %. The highest level it can reach is 1.3 volume %, and when that 
level is exceeded, the mission is called off. 

In the first scenario, we consider a case in which there is a leak of carbon 
dioxide from a carbon dioxide storage system. This means that CO 2 is constantly 
being added to the system and is likely to accumulate until some action is taken. 
Since this is unknown to the model, a good test of the filter would be to see how 
quickly this is detected. In this circumstance, it would be useful to monitor the 
levels and raise an alarm if the emergency levels are exceeded. Figure 5.5 shows 
how this situation leads to an emergency situation, when no mitigating actions are 
undertaken. The CO 2 concentration violates the SMAC at time, t = 320 s. 

In the second scenario, we simulate a fire in one section of the cabin, which 
is extinguished by a CO 2 extinguisher. (Halon cannot be used aboard the Space 
Station.) The CO 2 level consequently will immediately rise in the vicinity of the 
fire (both due to combustion product and due to the use of the extinguisher), and 
we wish to monitor how the level declines, and when the cabin becomes habitable 
again. This sort of a simulation would be invaluable in cases where there are multiple 
modules, and activity can be curtailed in the module under scrutiny until the levels 
are safe again. The release occurs at time t = 20 s, and continues for 40 s. Figure 
5.6 shows how the concentration of CO 2 changes with time for Sensor location #0. 
The SMAC is locally violated at t = 70 s, and the cabin is safe for habitation again 
at t = 200 s. The utility of the three-dimensional model lies in the fact that even 
local violations of safety standards can be detected, at both sensor and non-sensor 
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Figure 5.5: Average concentrations in the cabin in the presence of a continuous leak 
of CO 2 . The dotted-dashed line indicates the trror bounds outside which a fault is 
declared, and the dashed line is the SMAC levd for CO 2 . 
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locations even when the average concentration of the contaminant in the cabin lies 
below the SMAC. 

5.5.2 Other contaminants 

The monitoring of acute toxics would follow a slightly different procedure. 
For one, they are not normally present in the cabin atmosphere, so there are no sensor 
readings. Secondly, there are many specific toxics, each with different SMACS, and 
different sensors. Also, the tolerance for these substances will be tighter, so the 
procedure needs to be extra-sensitive. 

For this purpose, the best procedure would be to have a backup filter ready, 
and initialized, which can be activated as soon as any of the sensors register a reading 
for the toxic. The procedure can then access a central database for the SMAC for 
the substance, and begin to operate on the filter measurements. Since it is likely to 
be an unknown source, a diagnosis will have to be performed at the very beginning 
itself. A sample result is shown in Fig. 5.7. Sensor readings are identically zero, 
and the residual and the filter both show a zero reading. A source is introduced at 
time = 100 s , unknown to the filter. The filter responds almost immediately, and 
the fault is quite apparent within 10 s of the release. 

5.6 Summary and Conclusions 

Fault detection algorithms have been implemented using the error residu- 
als from the Implicit Kalman filter. The algorithm is able to detect and distinguish 
between sensor and process faults. The principle of Analytical redundancy using par- 
allel banks of filters is used to detect sensor faults, while process faults are detected 
when the residual estimation error from the filter exceeds pre-determined bounds. 
The filter is able to detect faults very quickly, which would be critical during space 
missions. 
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Contaminant tracking at Sensor location § 0 



Figure 5.6: Local concentration at a sample location in the cabin. The solid line 
indicated the true concentration at the Sensor location, and the dotted line, almost 
indistinguishable from the solid line is the filtered estimate of the concentration using 
the Implicit Kalman Filter 







Chapter 6 


Source identification - solving the inverse problem 


6.1 Introduction 

The final portion of this work is devoted to the inverse problem - that of 
identifying the unknown source that is causing a fault that has been detected by the 
detection algorithm. A rich body of literature exists in the realm of inverse problems 
(Alifanov, 1994; Kurpisz and Nowak, 1995), although much of the work has remained 
theoretical (Kirsch, 1996) with a few practical solutions. Part of the reason for this 
is the relative intractability of inverse problems, beyond simple cases with restrictive 
assumptions. 

One characteristic of an inverse problem is the unique manner in which 
data errors affect the error in the solution. The classic example given by Hadamard 
(Hadamard, 1923) was that of finding a soluticn u to the Laplace equation 

A„(x,,):=^> + ^l=0 i nM[0,co) (6.1) 

that satisfies the conditions 

u{x, 0) = f(x), — u(x, 0) = g(x), x € 0?, (6.2) 

where / and g are given functions. The unique solution for 


f(x) = 0 
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and 

g(x) = — sin(nx ) 

is given by 

u(x,y) = \sin{nx)sinh(ny),x G > 0. (6.3) 

n z 

With this solution, we therefore have 

sup{\f(x) + p(x)|} = 1 -4 0,n -> oo, (6.4) 

but then, for the error, 

sup\u{x, y)\ = \sinh(ny) — »• oo, n -> oo (6.5) 

n i 

for all y > 0. So, even though the error in the data tends to zero, the error 
in the solution it tends to infinity. 

Thus, a zero error in data tends to result in an infinite error in the solution 

it. 

6.2 Literature survey of solution methods 

Skliar (Skliar, 1996) used a one-shot optimization solution to estimate the 
location and capacity of a source, once it was detected. While this a relatively 
quick operation, and computationally non-intensive, it is prone to very high errors, 
especially in the presence of measurement noise. 

A study very similar in scope to ours was carried out (Richards et ah, 1997a; 
Richards et ah, 1997b) for application to a fire detection problem. They proposed 
a method for detecting, locating and sizing accidental fires in warehouses, based 
on the solution to an inverse heat transfer problem. They use a forward solution 
database, and minimize the least square error between estimated and measured times 
of activation of sensors that have been installed on the ceiling of the warehouses. On 
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closer examination, however, their problem turns out to be simpler since they assume 
a quiescent room with no air motion, and that the heat transfer occurs primarily 
through a buoyant plume of combustion gases that rises to the ceiling. Theirs is a 
one-shot solution mechanism. 

The estimation of multiple unknown sources is further complicated, and 
not many solutions exist. Cheng-Hung Huang and Jan-Yue Wu (Huang and Wu, 
1994) solved the two-dimensional inverse problem for two heat sources, but then, 
they assumed the boundary at which the sources were acting, therefore converting 
the problem to an inverse boundary problem. 

Mass transport inverse problems are common in geology where measure- 
ments at the earth’s surface are used to infer properties of processes occurring deep 
inside the earth. A general model has been developed (Talenti and Tonani, 1995) 
for gas-emitting geological systems, where the bulk gas velocity at the surface is 
used to locate the strength and location of tie gas source. Another application 
for inverse problems that dealt with mass transfer was developed by Australian 
researchers (Newsam and Enting, 1988; Enting and Newsam, 1990; Enting, 1993), 
who considered the problem of estimating surface sources of carbon dioxide and other 
trace contaminants from surface concentration data. They used a three-dimensional 
diffusion model for their transport process ana analytically solved the equation to 
account for the influence of various factors on the ability to invert measurement data 
to obtain source estimations. 

An elegant mathematical formulation for the determination of the source 
term was developed (Nanda and Das, 1996) for special cases of the heat conduction 
equation, which however assumes a specific mathematical form of the source function. 

6.2.1 Extended Implicit Kalman Filter 

Once a fault is detected, the next step is the identification of the source of 
the fault, namely its location and capacity. The 2-D model of Skliar used a one-shot 



105 


source identification, which is susceptible to high errors, especially in a system with 
so much noise (Skliar, 1996). Source identification problems of this type fall into the 
realm of inverse problems (Alifanov, 1994). An inverse problem , simply defined, is 
one where the cause is discovered from a known result. They arise in electrodynamics, 
geophysics, astrophysics and many other fields (Kurpisz and Nowak, 1995). Inverse 
problems axe ill-posed, which would mean that small data errors can lead to serious 
errors. Many fault diagnosis methods have been developed over the years, most of 
them applying to lumped systems. In general, they can be classified into pattern 
recognition (e.g. fault dictionaries), logic-based/information flow graphs (i.e., fault 
trees, signed directed graphs), and estimation/analytical redundancy methods. The 
reader is referred to one of many survey articles that address these issues (Basseville, 
1988; Frank, 1990; Gertler, 1991; Korbicz et al., 1991). While many tested techniques 
exist for lumped systems, distributed systems prove harder to solve because of the 
indirect relationship between the measurements and the model variables, and due to 
the large size of the model matrices. 

A multi step identification is proposed here, and the source identification 
will carried out over the time range between the time that a fault is suspected and 
the time that it is finally isolated. During the diagnostic process, the main filter will 
continue to run, but with a larger time step. At the end of the diagnostic process, 
the detected source term will be incorporated into the main mathematical model. 

The Implicit Kalman Filter, developed in Chapter 4, can be extended in 
order to estimate the unknown source. Extended Kalman filters are modifications 
to the Kalman filter (Halme and Seikainaho, 1986; Himmelblau, 1986) which can 
then be used to estimate both the state and parameters of systems. This is done 
through augmenting the state by adding the unknown source vector, f u to the state. 
Using a vector representation for the unknown source allows us to generalize the 
formulation to include unknown sources that are distributed spatially and multiple 


unknown sources. 
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The new augmented state is now 


Q 


a 

m 



( 6 . 6 ) 


where Q m , the augmented state is given by Equation 4.5, and f u , the unknown 
source vector of dimension n is given by the equation 


fu = 


0 

0 

0 

0 


fy-i 

0 

0 


(6.7) 


Here, f ui represents the capacity of a single unknown source. The vector, f u 
can be modified to handle a single source distributed spatially, or multiple sources. 

In addition, we assume that the source term is relatively unch an ging, and 
that it satisfies the equation 


-^ = 0 + w !lm (6.8) 

where w um is the uncertainty associated with he unknown source. 

Integrating and discretizing the equation, we obtain the following equation 
for the description of the new state. 

fum+l = fum + kmW um (6.9) 


where k m is a matrix reflecting the integration ,ime step. The noise in the unknown 
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source is characterized by its covariance, S m , which is a user-defined constant that 
depends on how much the unknown source capacity can vary. High values for S m 
reflect a high uncertainty in the unknown source capacity, and will result in a filter 
that responds quickly to the residual errors, whereas low values for S m reflect a fairly 
constant unknown source capacity and will result in a filter that responds slower to 
the residual errors, but one that will provide a smoother solution for the unknown 
source estimation. 

The state equations, once the unknown source is included now read 


Ai a Qm+i = A 2 a Q^ + 

[ 

1 


r -I 


W m 

+ 

Cm km 



. 0 




Wum 


( 6 . 10 ) 
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(6.11) 


0 0 (A x /2 + A y + A z — r/At) W 


A 2 “ = 


0 0 
0 0 
0 0 


Ay/2 

— A z /2 

0 


0 

0 

w 


(6.12) 


W is an n x n diagonal matrix, where n is the number of grid points used for the 



108 


discretization, and has the structure, 


0 


W = 


0 


0 


0 


1 

0 o 


(6.13) 


0 


The matrix is extremely sparse, with all zeros except for source location, where there 
is a unit term. This unit term picks off the source capacity scalar term on expansion 
to yield Eq. 6.9. For multiple sources or single sources that are distributed spatially, 
there will be as many non-zero terms as there are sources, placed at the appropriate 
grid points. 

The measurement matrix remains unchanged since we have no measure- 
ments that have a direct bearing on the unknown source, but the matrix has to be 
augmented. 


z m+ 1 


Hi(m + 1) 0 


Qm+l v m+l 


(6.14) 


The filter changes to account for the new state and the estimation algorithm 
consists of an Implicit and an explicit part, the Implicit part used in the estimation 
of the state, and a parallel explicit portion usee in the estimation of the strength of 
the unknown source. The gain terms serve to n ultiply the residual errors generated 
in order to update the state estimation for the lext time step. 
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Filter (Implicit) : 


(-** - — )q' 

1 2 A t ,H 
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(6.15) 


Filter(Explicit): 


fum+l — ^um "k ^4[z Hi y m-{- 1 [m] 


(6.16) 


The predicted estimated value for the auxiliary variable y a is given by the 
following equation 


y m+l|r 


a 2 


fm + fum 

w 

4m|m "f" 

0 


(6.17) 


where y a is the new augmented estimate and is a column vector of size 4n. 

The augmented Implicit Kalman Gain matrix is a 4n by m matrix, with 
four partitions (one for each co-ordinate direction and one for the unknown source) 
and has the structure 


L 


a 

m+1 — 


T T T t T 

IjI Li 2 



and is obtained from the equation 


(6.18) 
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(6.19) 
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where 

P m+i|m a = (6.20) 

' A5 3 P’ |m Af T + A5>P’ |m Al» T A^P’ |m Al> T P“ |m + CQS* 

+ p m|m + CQC T 

AfP’| m Aj 3T AfP’ |m Af T AfP’ |m A|» T 0 

A 2 3P mlm A 2 3T A i 3 P q | Af T Ai 3 P q , Ao 3T 0 

P m|m + S Q CT 0 o P m|m + S QS T 

P m+i| m a is the augmented predicted error covariance matrix, and is of 
dimension 4n by 4n. P^| m is the uncertainty matrix associated with the unknown 
source, and in general is an n by n diagonal matrix, with the uncertainty due to each 
unknown source location as the diagonal elements, though it reduces to a sparse 
matrix with just one diagonal element, for a fault caused by an emission at one 
spatial location. P^,| m is the covariance of the estimated error, and defined by Eq. 
4.21. C and Q represent the stochastic model disturbance transition matrix and the 
model covariance matrix, respectively. The superscripts in the A 2 matrix refer to 
the appropriate partitions of the A 2 “ matrix defined earlier. 

The gain matrix partitions, L 2 and L 3 remain unchanged from the simple 
implicit filter, and only the Li partition changes, along with the introduction of the 
new gain partition for the unknown source, L 4 The gain matrix partition, Li now 
becomes a function of the uncertainty term of the unknown source, P u , which is 

1 m| m ’ 

then used in the correction term for the first Implicit Kalman Filter equation in Eq. 
4.13. 

Example: 

For a system with six sensors and one unknown source, which is the most 
commonly expected fault, the gain matrix partition for the unknown source term 
has the structure 



Ill 



where i refers to the location of the unknown source in the cabin, i 4 ^ is a row vector 
of dimension m, and is evaluated by simplifying Eq. 6.19 for one unknown source. 
In general, L 4 has as many non-zero rows as there are unknown sources. 



The lower case symbols, p, s, q , and c represent the scalar terms for the uncertainty 
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associated with the unknown source, the transformation factor representing the inte- 
gration term in the unknown source terms, and the model uncertainty for that grid 
location. 

The measurement term for the grid point where the unknown source is 
located for sensor s. is The computation therefore yields I4 ( , which is a row 

vector of dimension 6, since we have six sensors in this case. Each of these six gain 
terms multiplies the appropriate sensor residual from z - Hiy , which is then used 
to update our estimate for the unknown source location via Eq. 6.16. 

The gain term, in turn is used to propagate the uncertainty for the next 
time step, which is governed by the equation 

P y ™+I| m+ 1 = [I - I4+iHt]P" +l|m (6.23) 

and 

P y m+l|m+l = > (6.24) 

6.2.2 Initial guesses through sensitivity matrices 

The Extended Implicit filter, like all Kalman filters requires an initial guess 
for fuT7ii which is not a trivial problem. While the filter performs well in filtering 
out noise, and adjusting for model errors, its performance in estimating an unknown 
source depends crucially on a good initial guess for the location of the unknown 
source emission. 

For the purpose of this derivation, we assume that there is only a single 
localized source term that is causing the fault. The solution can be extended to 
multiple and distributed sources, but it would be more complicated since the single 
source will be broken down into a combination of linearly independent single sources. 

In this work, we use pre-calculated sensitivity coefficients for the purpose. 
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Sensitivity coefficients have been widely used in estimating solutions to inverse prob- 
lems, especially in areas of heat conduction, geology, and tomography (Alifanov, 
1994). Some detailed analysis of sensitivity coefficients and their properties and use 
are also available (Beck et ah, 1985; Kurpisz and Nowak, 1995). The sensitivity 
coefficient can be defined as 


Z = 


dq 

df u 


(6.25) 


where f u represents the source term. The sensitivity coefficient term represents the 
sensitivity of the concentration q at each mesh point with respect to the release of 
a source at every mesh point. Z is calculated by solving the basic model equation, 
Eq. (3.1). While Z is defined over the entire domain, the only points of interest 
will be the sensor locations since those are only points about which we have direct 
information about the concentration that we can use in the event of a fault. Z is 
therefore partitioned into a usable and non-usable part, 


Z = 


Z 


sensors 


z 


non— sensors 


(6.26) 


This pre-calculated sensitivity matrix can be computed for different times, 
in order to provide a window of time over which the fault diagnosis can be conducted, 
depending on how soon the fault is detected. 

The first step, therefore, is to calculate Z. In order to do this , we multiply 
Eq. 3.1 throughout by £■, which leads to the equation 


d 

df u 


dq 

dt 


+ u • Vq = Dm V 2 q + f u 


(6.27) 


Since f u is independent of the co-ordinate axes, we can rewrite Eq. 6.27 as 


d dq 
dt df u 


+ u • V 


«)- 


(6.28) 


Replacing Z for JyjL from the definition for the sensitivity coefficient, we 
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get the direct well-posed problem equation 

dZ , 

— + u • VZ = D m V 2 Z + 1 (6.29) 

subject to the appropriate boundary and initial conditions. The identity matrix, I 
essentially is a unit source that is sequentially placed at every grid location in the 
cabin in order to measure its effect on all the other grid locations. The boundary 
conditions and the initial conditions will also have to be divided throughout by ^ 
to obtain the initial and boundary conditions for the sensitivity problem. The Z 
matrix is then partitioned to obtain Z sensors . Since the structure of the equation 
is unchanged from the original model equation, the same algorithms and numerical 
techniques can be used in computing the sensitivity matrix. For each time step, the 
computation of the sensitivity matrix takes about 180 CPU minutes on the DEC- 
Alpha 500 AU, and therefore, the computations have to be performed off- line and 
before the time that the diagnosis can take place. 

The solution to Eq. 6.29, Z (t), is a function of time, where the time refers 
to the time elapsed since the fault occurred. For the purpose of using the Z sensors 
partition in calculating an initial guess for the capacity and location of an unknown 
source, we use 


Z critical — ’^‘sensort |t=td e t„t (6.30) 

where tdetect refers to the time when the fault is detected. This yields a Z critlca i matrix 
of dimensions n x m, where n is the number of grid points, and m is the number of 
sensors. In other words, the Z critical matrix contains the response observed at each 
sensor location for a unit perturbation at every location in the cabin. In the event of 
a fault, this response is used in conjunction with the observed measurement response 
in order to estimate the perturbation that cam ed the fault. 

The solution proposed here consists of two parts; a first off-line part that in- 
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volved pre-calculating sensitivity matrices, and which is a computationally intensive 
process; and a final real-time computationally non-intensive portion that actually 
computes a first guess for the location and the capacity of the unknown source emis- 
sion. 

The strategy of the method is as follows: 

• The ill-posed problem is made well-posed by making assumptions about the 
problem in the areas of ill-posedness, in this case the unknown source term 

• The well-posed direct problem is solved for this assumed value(s). 

• The measured quantities are noted for the ill-posed problem using the sensor 
system 

• The calculated values for the assumed problem are compared with the mea- 
sured values from the sensor system, and the assumed input data are modi- 
fied to ensure a matching of these two quantities. Although this method is 
likely to yield a correct solution, the nature of an ill-posed problem could 
mean that there are multiple solutions, with the method being able to iden- 
tify only one of them. 

Once the sensitivity matrix has been calculated, we now demonstrate how 
it is used in the fault diagnosis process. 

In the event of a fault being diagnosed, an error m-tuple is generated using 
the prediction errors from the Implicit Kalman Filter. 

ci 
C2 
e3 

e m— 1 
Cm 


e = 


(6.31) 
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where e\, e 2 , e m _! and e m represent the prediction errors for Sensors l,2,m-l and m, 
and represent the error between the projected estimate and the measurement signal, 

e i = z i~ J/im+lIm (6.32) 

Next, the capacity of the unknown source is calculated, sequentially assum- 
ing that the emission occurred at each location in the cabin, 

2 - 

Cap { = — (6.33) 

e ? 

This calculation results in n projected capacity m-tuples, with the assump- 
tion that the source is in each of n locations, where n refers to the number of grid 
points being used. The algorithm being used for this estimation (Fig. 6.1) is given 
below. 

(1) Calculate capacities using Eq. 6.33. 

(2) Scale capacities within each m-tuple using the maximum and minimum ca- 
pacities within the m-tuple. 

(3) Compute the standard deviation of the calculated capacities for each m- 
tuple, for both of the scaled versions. 

(4) Pick the point with the least standard deviation. If the same point is picked 
using both the minimum and maximum scalings, that point is likely to be 
the source location. 

(5) If minimum and maximum scalings yi Bid different locations, use the loca- 
tion which yields the lower standard deviation as the starting guess for the 
location of the unknown source. 

(6) The initial guess for the capacity is calculated from Step 1. 
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Figure 6.1: Sensitivity analysis algorithm that yields a guess for the location and 
capacity of an unknown source 
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Table 6.1: Me asurements after 1 time step 


Source strength 

Measurement 

0 

0.3000001 

1 

0.3005528 

50 

0.3276336 

100 

0.3552671 


The final estimate is based on the principle that for the correct assumption 
of the source, the m-tuple of the projected capacities would show the minimum 
standard deviation, and would match the previously generated m-tuple of sensitivity 
coefficients for the location, to a multiplicative real constant. 

This is an area of work that is amenable to treatment using Artificial In- 
telligence or Knowledge Based systems, where the algorithm is trained to pick up 
patterns and suggest intelligent solutions based on learned past experience. 

6.2.3 Sensitivity experiments 

Tables 6.1 through 6.5 are the concentration measurements at Sensor lo- 
cation #6 (8,14,18) for a unit source release at location (8,15,12). This location 
happens to be at a sensor location closest to this source, and consequently has a 
high response. The data are shown to indicate how the sensitivity indices helps 
scaling. 

The experiment was repeated for a se isor location farther from the source 
emission site, location #3 at (3,5,7). As expected, the sensitivity coefficients are 
smaller in magnitude. Tables 6.6 through 6.10 show the concentration measurements 


Table 6.2: Measurement at Sensor location #6 after 2 time steps 


Source strength Measurement 


0 

0.3000004 

1 

0.3013975 

50 

0.3698579 

51 

0.3712551 

100 

0.4397156 
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Table 6.3: Measurement a t Sensor location #6 after 10 time steps 


Source strength 

Measurement 

0 

0.3000708 

1 

0.3147669 

50 

1.034878 

51 

1.049574 

100 

1.769685 


Table 6.4: Measurement a t Sensor location #6 after 100 time steps 


Source strength 

Measurement 

0 

0.2533103 

1 

0.2858888 

50 

1.882236 

51 

1.914815 

100 

3.511162 


Table 6.5: Measurement a t Sensor location #6 after 300 time steps 


Source strength 

Measurement 

0 

0.121 

1 

0.1617357 

50 

2.115451 

51 

2.155323 

100 

4.109037 

1000 

39.99361 
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Table 6.6: Measu rements at #3 after 1 time step 


Source strength 

Measurement 

0 

0.2999996 

1 

0.2999996 

50 

0.2999996 

51 

0.2999996 

100 

0.2999996 I 


at (3,5,7) for different times. 

One interesting result is that the source is not detectable by this sensor 
during the first two time steps, and even at 10 times steps, the effect is very slight. 

The distance between the source emission location and the sensor location 
is efficiently captured in the sensitivity coefficient, which is computed earlier. The 
use of the sensitivity matrix in estimating the detectability of a fault is discussed 
later. 

6.2.4 Sensitivity analysis-Results 

In a real setting, the sensitivity analysis would be tested by introducing 
a fault, taking sensor readings, and checking to see if the sensitivity analysis will 
estimate the source correctly from the readings. In the absence of an experimental 
setting, we introduce the unknown fault, run the model, and use the model calculated 
concentrations at sensor locations as sensor readings. In order to simulate a real 
physical setting, we also add a Gaussian nois ; to the sensor readings to generate 
pseudo-experimental values for testing the Sensitivity analysis. 

For our test case, we first pre-calculate the sensitivity matrix for the sensor 


Table 6.7: Measurement at Sensor lo< ation #3 after 2 time steps 


Source strength 

Measurement 


0.2999778 

i 

0.2999778 


0.2999778 

51 

0.2999778 

100 

0.2999778 







Table 6.8: Measurement a t Sensor location #3 after 10 time steps 


Source strength 

Measurement 

0 

0.2826567 

1 

0.2826571 

50 

0.2826755 

51 

0.2826759 

100 

0.2826943 


Table 6.9: Measurement a t Sensor location #3 after 100 time steps 


Source strength 

Measurement 

0 

0.1181842 

1 

0.1196233 

50 

0.1901395 

51 

0.1915785 

100 

0.2620947 


Table 6.10; Measurement a t Sensor location #3 after 300 time steps 


Source strength 

Measurement 

0 

0.044009563 

1 

0.048030302 

50 

0.2450462 

51 

0.2490669 

100 

0.4460829 
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Table 6.11: Sensor location and the associated measurement noise 


Sensor 

Co-ordinates 

R 

1 

10 10 10 

0.0001 

2 

8 2 23 

0.001 

3 

3 5 7 

0.0001 

4 

6 6 6 

0.0023 

5 

4 11 22 

0.019 

6 

8 14 18 

0.0019 


array shown in Table 6.11. The base-line concentrations are taken to be the steady 
state concentration field for an initial CO 2 concentration of 0.6 volume %. 

A source of strength 500 mg/m 3 is then introduced into the model, acting 
at location (8,6,5), and the pseudo- measurements that are generated from the model 
are used for the sensitivity analysis. A pseudo-random vector with six elements with 
a bound of 5 % of the measurement was generated and added to each measurement 
in order to mimic a sensor system with 5 % measurement noise. 

The sensitivity analysis reports location (8,6,6) as the one with the least 
standard deviation (ct = 23.67), and the capacities as calculated by each of the 
sensors are reported in Table 6.12, the mean of which will serve as our initial guess 
for the Extended Implicit Kalman Filter. 

While the sensitivity analysis’ initial estimate is quite satisfactory for the 
purpose of the filter, we also note that the next highest standard deviation (a = 
200.86) was reported for the correct location (8,6,5), whose guess capacities are 
reported in Table 6.13. One observes that Sensors 1-4 have excellent estimates for 
the capacity, but the correct location is not reported by the sensitivity analysis 
since Sensors 5 and 6 report capacities that increase the standard deviation of that 
particular m-tuple. 

Eq. 6.29 yields Z sensors that can be evaluated for different times. This can 
provide a window of time over which the fault diagnosis can be conducted, depending 
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Table 6.12: Estimated capacities for assumed source at (8,6,6) 


Sensor 

Capacities (mg/m 3 ) 

i 

158.2 

2 

140.7 

3 

154.7 

4 

158.6 

5 

103.8 

6 

115.3 


Table 6.13: Estimated capacities for assumed source at (8,6,5) 


Sensor 

Capacities (mg/m 3 ) 

i 

501.9 

2 

512.5 

3 

512.5 

4 

512.5 

5 

105.5 

6 

137.4 
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on how soon the fault is detected. In our examples, we have only used Z critical, and 
this could cause large errors if the fault were not accurately detected. For cases where 
one can expect a substantial delay in the detection of a fault, it would be prudent to 
minimize the standard deviation of the capacities, evaluated over a period of time, 
and using say, Zj = i, Z^-iq, and Z$ = 5 q, which would provide added robustness to the 
analysis. 

6.2.5 Iteration procedure for Extended Implicit Kalman filter 

As we observed in the previous section, the sensitivity analysis can some- 
times pick a point in the vicinity of the actual source, instead of the actual point 
itself. The presence of turbulence in the cabin air, and general measurement noise 
can increase the distance of the actual solution from the first guess. In order to refine 
our solution, we use the property of the Implicit Kalman Filter being an optimal 
estimator and run the filter for the initial guess and determine the squared prediction 
error. The final solution is reached when the point under consideration has the least 
squared error when compared with its six nearest neighbors. 

The algorithm for this location search is shown in Fig. 6.2 and has the 
following steps. 

(1) Choose one of B’s spatial neighboring points (grope), which has not been 
previously visited as the assumed location for the unknown source. 

(2) For this location, and the original estim ite for the capacity, run the Extended 
Filter and obtain a new squared predic:ion error, and an estimated capacity 
for that location. 

(3) Repeat Steps 1 and 2 until all neighbor; not previously considered have been 
covered. 

(4) Compare the Squared prediction error, and for the next approximation, 
choose the point with the least SPE. 
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Figure 6.2: Extended Kalman Filter procedure for estimating unknown source 
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Table 6. 14: SPE for neighbor points for a converged solution 


Co-ordinates of neighbor points 

SPE 

7,6,9 

64.557 

7,5,9 

293.19 

7,7,9 

597.72 

6,6,9 

9761.09 

8,6,9 

40655.4 

7,6,8 

1581.80 

7,6,10 

3070.07 


(5) Repeat Steps 1 through 4 until a minimum is obtained. 

6.3 Results 

In the first test case, a source of strength 500 mg/m 3 is introduced at 

(7.6.9) , with an initial guess of 400 mg/m 3 for the capacity at location (7,6,9) from 
the sensitivity analysis. The Extended Filter converges rapidly to the final solution 
(Fig. 6.4), and the SPE is only 64.557. While the solution appears to be acceptable 
on inspection, our algorithm requires that the predicted errors be examined for its 
neighbors. Accordingly, we run the filter for all the neighbors to (7,6,9), and confirm 
from Table 6.14 that the estimate is indeed the best, given the measurement data. 

In the next test case, a source of strength 1500 mg/m 3 was applied at loca- 
tion (11.11,9), and the measurements generated were used as inputs to the sensitivity 
analysis algorithm. For this test case, the sensitivity analysis points to a source lo- 
cation at (12,10,8), and a capacity of 1150 mg/m 3 . The filter does not converge for 
this location, and so we try the neighbors. 

Based on the squared error results in Table 6.15, we shift our focus to 

(12.10.9) since it has least error and rerun the filter for its neighbors. 

In two iterations, we are able to pinpoint the correct location and source 
capacity in the presence of noise. Figure 6.8 shows how the squared prediction error 
converges in three iterations to its final minimum value. 
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Figure 6.3: Comparing the squared estimation error across nearest-neighbors 


Table 6.15: SPE for neighbor points for a converged solution-iteration 1 


Co-ordinates of neighbor points 

SPE 

13 10 8 

10545.7 

11 10 8 

7719.28 

12 11 8 

1.384 X 10 6 

12 9 8 

554869.0 

12 10 7 

552.204 

12 10 9 

157.859 
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Figure 6.4: The estimation of a constant source (unknown to the model) of 500 
mg/m 3 by the Extended Implicit Kalman Filte *, with an initial guess of 400 mg/m 3 . 


Table 6. 16: SPE for neighbor points for i. converged solution-it eration 2 


Co-ordinates of neighbor points 

sum-squared error 

12 11 9 

72.6206 

12 9 9 

429.040 

12 10 8 

233.676 

12 10 10 

195.053 

11 10 9 

10458.9 

13 10 9 

10158.0 
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timesteps 


Figure 6.5: Unknown source capacity estimation - The first of a series of three 
iterations needed for the Extended IKF to converge to its final solution. Note that 
because of the wrong guess for the unknown source, the filter (wavy line) under- 
predicts the actual strength of the source (straight line). 
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timeste is 


Figure 6.6: Unknown source capacity estimation- The second in a series of three 
iterations needed for the Extended IKF to converge to its final solution. Note that 
the estimation is better in this case than it is for Iteration 1, but the filter (wavy 
line) still slightly under-predicts the actual str<ngth of the source (straight line). 
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Figure 6.7: Unknown source capacity estimation- The third and final Extended 
IKF iteration results showing the estimated (wavy line) and correct unknown source 
(straight line) capacities for the correct guess for the unknown source location. 
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If for the same case, i.e., for a source of strength 1500 mg/m 3 introduced 
at (12,11,9), an initial guess of 1400 mg/m 3 for the capacity from the sensitivity 
analysis is used, the Extended Filter converges rapidly to the final solution (Fig. 
6.9), and the SPE is only 22.95. The quality of the initial guess, therefore, is very 
crucial to the performance of the filter in terms of a quick convergence to the final 
solution. 

The Extended Filter is only as good as the model and the sensor system 
are. If all the sensors are crowded in one area of the cabin, faults in other regions 
can go undetected. 

As an example, consider a fault at location (13,10,20). A constant source 
of strength 500 mg/m 3 is applied over a time of 300 time steps, and the sensitivity- 
analysis was used to estimate the location and capacity for the unknown source. 
The guess location was a significant distance away from the correct location. The 
filter was then run on these measurements, assuming the correct location as the 
guess location with an initial guess of 400 mg/m 3 . Fig. 6.10 shows the tracking 
for this case, and it indicates that the response of the filter has been too slow, 
and even though the filter is starting to approach the final solution, the solution 
does not converge within the desired duration of time. It is therefore crucial to be 
aware of ‘dead-zones’ which are outside the domain of observability and consequently, 
detectability. 

6.3.1 Varying functions 

The Extended Filter can be applied even to varying sources. If the sources 
are varying very quickly, then our assumption is violated, and therefore the solution is 
likely to be oscillatory and posses a high residual error. Figure 6.1 1 is an example of a 
varying source that was tracked by the filter. The filter exhibits oscillatory behaviour 
because the unknown source function is changing rapidly and because it operates 
with no knowledge about the nature of this unknown source function. 
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Figure 6.9: Unknown source capacity estimation- A very good first guess can lead 
to the Extended IKF rapidly converging to the correct solution. 
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Figure 6.10: Unknown source capacity estimation when the fault is outside the active 
observability range 
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Figure 6.11: Varying source capacity estimation. The Extended Filter is able to 
diagnose rapidly varying functions, though wit l some overshoot and undershoot. 
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In another test case, a source function that varies more slowly is used. 
Figure 6.12 shows the input and estimated source functions. The overshoot and un- 
dershoot is unavoidable since the filter has no a priori knowledge about the function. 

If there were any information about the nature of the source function; for 
example, if a cylinder had sprung a leak, and the escape of contaminant was governed 
by a certain equation, Eq. 6.9 could be modified to reflect that information, and the 
filter would be modified accordingly. In the absence of any information, the current 
assumption that the source remains constant or slowly varies is found to produce 
satisfactory results. 

6.4 Sensitivity analysis-maps 

The calculation of the sensitivity coefficients, although time consuming is 
crucial to obtaining a good initial guess for the location and capacity of the unknown 
source. In addition, it also serves as a useful tool in determining the observability, 
and consequently the fault detectability, of the system. 

The system under study presents some unique features with respect to its 
observability and controllability. The usual methods of determining observability 
and controllability (Ramirez, 1994) fail here. In the strict sense, if the state of the 
system is defined to be the concentrations at each mesh point, and if the controls 
involve varying the inlet and outlet velocities and their concentrations, this system 
is neither controllable nor observable. This is true for many distributed-parameter 
systems, especially of the reaction-diffusion parabolic type because of the way a 
perturbation propagates through the system, and for systems where the noise can 
be substantial. The sensitivity matrix contains useful information about the speed 
with which a perturbation travels through the system, and can be used to plot zones 
of observability. 

Visualization of the zones would be critical in fault detection and diagno- 
sis, especially for real-time applications, so that the operator can easily ascertain the 
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certainty of the diagnosis. The zone is visualized through contours of equal sensi- 
tivity coefficients. Grid points that are along a constant sensitivity coefficient have 
equal observability and detectability. The sensor uncertainties can then be used to 
determine the lower limit of the sensitivity coefficient that can be detected. Grid 
points that have Z values below the threshold cannot be detected, and lie outside 
the contoured zones. 

Figures (6.13-6.14) show these zones for two sample slices across the cabin. 
The observable zones lie within the region bounded by the contours. The overlap 
zone of the six sensors is clearly visible. Experiments also show that faults occurring 
outside the zone are difficult or impossible to diagnose accurately, though an accurate 
estimation of unknown sources in this zone is sometimes possible when the sensor 
uncertainties are very small compared to the model uncertainty (accurate sensors). 
The observability contours can also serve as an excellent tool for sensor selection and 
placement because the observability contours for different sensor configurations can 
be visualized in order to achieve maximum coverage of the module. 

6.5 Summary and Conclusions 

In this Chapter, a method for tackling a specific kind of inverse problem 
was developed. Inverse problems are usually ill-posed, posses no unique solutions, 
and small data errors can cause very large errors in the estimated solution to the 
problem. Most inverse problems tackled in the literature are boundary-value inverse 
problems, and sometimes assume a certain mathematical functional form for the 
function that is to be estimated. The inverse problem that arises out of the attempt 
to estimate the capacity and location of an unknown source provides the added 
complication that it is not a boundary-value problem, and that fact that there is 
little or no information about the nature of the source. Our solution technique 
consists of a combination of two commonly used techniques, sensitivity analysis, and 
Kalman filtering, which provides a very effective tool for estimating the location and 
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Figure 6.13: Observability and detectability coitours for slice at grid height 15 (lm 
above the floor of the cabin). The region within the contours is the observable zone. 
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Figure 6.14: Observability and detectability contours for slice at grid height 25 (1.7m 
above the floor of the cabin). The observable zone lies within the contours. 
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capacity of the unknown sources. The two methods individually would be ineffective 
since sensitivity analysis does not respond very well to noisy data, and Kalman 
filtering is too computationally intensive to use without a good first guess for the 
source location. The combination of the two methods should be further investigated 
for its properties, and could possibly be extended to many applications that involve 
distributed parameter systems. 



Chapter 7 


Concluding remarks 


7.1 Original contributions 

The focus of this work was to come up with implementable algorithms for 
use in Space missions. This research has shown that advanced control techniques can 
be readily modified for use for this real-time application. Specifically, we have shown 
that the detailed modeling of the dispersion of air contaminants, their monitoring and 
detection can be accomplished with reasonable computational power. Every attempt 
has been made to tailor the work to NASA specifications (Technology Development 
Requirements, 1996). 

Among the original contributions of this work are: 

• Proving the inadequacy of lumped systems for air contaminant monitoring, 
and the limitations of two-dimensional models. 

• Using CFD techniques for modeling the flow aboard the Space Station, hith- 
erto not attempted 

• Development of the Extended Implicit Filter, and using it in combination 
with Sensitivity analysis for joint state and parameter estimation. 

• The use of sensitivity maps in identifying dead zones, and observable/detectable 
zones aboard the Space Station. 
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• Overall, providing a framework for accurate and efficient air contaminant 
monitoring using algorithms and methods suited for real-time implementa- 
tion, and including sensor placement decisions. 

7.2 Possible applications 

This research can be used as a starting point for implementing future gen- 
eration air contaminant monitoring systems for the Space Station and for long range 
Space missions. Some of the specific tasks that this research can be used for include: 

• Design of habitation geometry and forced flow sources and sinks. This re- 
search can be used to identify stagnation zones in the cabin, and to design 
better ventilation systems to maximize crew safety and comfort. 

• Sensor selection and placement issues- The sensitivity coefficients have been 
shown to be an effective measure of the observability and detectability of 
faults aboard the spacecraft. The visualization of the sensitivity coefficient 
provides an easy tool for facilitating sensor selection in order to ensure that 
all parts of the cabin are adequately covered by the measurement system. 

• Analysis of proposed responses to emegencies aboard the spacecraft. The 
algorithms developed here can be used to simulate possible accidents, and 
study the effect of the remedial measures on the spacecraft atmosphere. For 
example, one could study how long it would take for a toxic leak in one 
module to spread to another module, and how long it would take for a 
module to be habitable after a release ( ccurred and the removal devices are 
activated. In addition, proposed remediation measures for many different 
kinds of accidents can be analyzed to ascertain their relative merits and 
demerits. 

• Analysis and design of fire extinguishment systems. The effect of fires, the 
transport of combustion products, and t he air flow patterns in the likelihood 
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of the air supply being cut off can be studied for a variety of conditions. 

• Design of automated and semi-automated control systems, which use the 
measurement system and model inputs to provide for operations that do not 
need human intervention. 

7.3 Directions for future research 

Computational Fluid Dynamics (CFD) has not been adequately utilized in 
Space station environmental control applications. CFD could prove to be a very 
powerful tool in the design of Space station modules for crewed missions. A detailed 
study of these flows is therefore in order. The next step needs to be an experimental 
validation of our results under both atmospheric and micro-gravity conditions. The 
validation needs to occur in two specific areas-one is to study room air flows, using 
the cabin if possible, and observing how much noise is usually present in the mea- 
surements, and in noting the performance of the filter itself. Tuning the gains via 
sensor and model uncertainties is another area of work that needs scrutiny. Simula- 
tions have shown that a proper choice of these parameters will affect filter gains, and 
the consequent performance of the filter. The Implicit Kalman Filter and Extended 
Filter could both be optimized to improve performance. Alternative formulations of 
the filter which do not require the inversion of matrices would reduce the number of 
operations needed per time step. 
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Nomenclature 

= State Transition Matrix (Left Hand Side) 

= State Transition Matrix (Right Hand Side) 

= Discrete representation of the spatial operator for the State 
Transition Matrices in each of the co-ordinate directions 
= Zero-mean white Gaussian processes 
= Dimensionless empirical constants used in the 
k — e turbulence model 

= Stochastic model disturbance transition matrix 
= Discrete analog of the mass or eddy diffusivity ( m 2 s~ l ) 

= Mass Diffusivity (m 2 s~ x I 
= Eddy Diffusivity (m 2 s -1 ) 

Error m-tuple used in sensitivity analysis 
= Contaminant source capacity (kgm~ 3 s^ 1 ) 

— Unknown contaminant source capacity (Jt^m” 3 s _1 ) 

= Body force per unit masi acting on fluid N/m 
= Elements in the measurement matrix, H 
= Measurement matrix 
= Identity matrix 

= Kinetic energy of turbulence (used in the k — e model) 

= Constant of integration lor the unknown source capacity 
= Loss function 
= Implicit Kalman Gain 
= Pressure of fluid Nm ~ 2 
= Predicted error covariance matrix 
= Estimation error covariance matrix 
= Covariance matrix of the state 
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Subscripts 
n, p, r 
Superscripts 
a 

E 


= Concentration of Contaminant (fcpm“ 3 ) (volume %) 

= Discrete analog of contaminant concentration 
(fcpm -3 ) or (volume %) 

= Volumetric heat flux addition/removal to the fluid Jm' 3 s _1 
= Covariance matrix (diagonal matrix) of the model noise 
= Concentration vector 
= Row unity matrix 

= Covariance matrix (diagonal matrix) of the measurement noise 
= Covariance matrix for the unknown source capacity noise 
= time (s) 

= Fluid temperature (K) 

— Velocity vector 

— Velocity components in the co-ordinate directions (ms~ l ) 

— Stochastic model noise 

= Location(s) of the unknown source(s) used in fault diagnosis 
= State vector 
= Estimate of State 
= Estimation Error 

= Co-ordinate directions or positions (m) 

= Measurement signal 
= Sensitivity matrix 

= Mesh indices in the coordinate directions 

= Augmented quantities, augmented to include unknown 
source(s) in the model 
= East 


W 


= West 
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N = North 

S = South 

U = Up 

D = Down 

Greek Symbols 

Ox, «z = State transition matrices used in the double filter 

^7 = discretization steps in the spatial directions (ADI method) 

e = Dissipation rate of the turbulence (used in the k — e model) 

A = Second viscosity coefficient (m 2 s -1 ) 

fi = Molecular Viscosity coefficient (m 2 s _1 ) 

Ht = Eddy Viscosity coefficient (m 2 s _1 ) 

&k = Error Functional 

p = Density of fluid ( kgm~ 3 ) 

°k ■ a i — dimensionless empirical constants used in the turbulence model 



